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ABSTRACT 

We present a self-consistent model of the spectral energy distributions (SEDs) of spiral galaxies from the ultraviolet (UV) to the 
mid-infrared (MIR)/far-infrared (FIR)/submillimeter (submm) based on a full radiative transfer calculation of the propagation of 
starlight in galaxy disks. This model predicts not only the total integrated energy absorbed in the UV/optical and re-emitted in the 
infrared/submm, but also the colours of the dust emission based on an explicit calculation of the strength and colour of the UV/optical 
radiation fields heating the dust, and incorporating a full calculation of the stochastic heating of small dust grains and PAH molecules. 
The geometry of the translucent components of the model is empirically constrained using the results from the radiation transfer 
analysis of Xilouris et al. on spirals in the middle range of the Hubble sequence, while the geometry of the optically thick components 
is constrained from physical considerations with a posteriori checks of the model predictions with observational data. Following the 
observational constraints, the model has both a distribution of diffuse dust associated with the old and young disk stellar populations 
as well as a clumpy component arising from dust in the parent molecular clouds in star forming regions. In accordance with the 
fragmented nature of dense molecular gas in typical star-forming regions, UV light from massive stars is allowed to either freely 
stream away into the diffuse medium in some fraction of directions or be geometrically blocked and locally absorbed in clumps. 
These geometrical constraints enable the dust emission to be predicted in terms of a minimum set of free parameters: the central 
face-on dust opacity in the B-band t b , a dumpiness factor F for the star-forming regions, the star-formation rate S FR, the normalised 
luminosity of the old stellar population old and the bulge-to-disk ratio B/D. We show that these parameters are almost orthogonal in 
their predicted effect on the colours of the dust/PAH emission. In most practical applications B/D will actually not be a free parameter 
but (together with the angular size 9 gaj and inclination i of the disk) act as a constraint derived from morphological decomposition 
of higher resolution optical images. This also extends the range of applicability of the model along the Hubble sequence. We further 
show that the dependence of the dust emission SED on the colour of the stellar photon field depends primarily on the ratio between 
the luminosities of the young and old stellar populations (as specified by the parameters S FR and old) rather than on the detailed 
colour of the emissions from either of these populations. The model is thereby independent of a priori assumptions of the detailed 
mathematical form of the dependence of SFR on time, allowing UV/optical SEDs to be dereddened without recourse to population 
synthesis models. 

Utilising these findings, we show how the predictive power of radiative transfer calculations can be combined with measurements of 
Bgai, i and B/D from optical images to self-consistently fit UV/optical-MIR/FIR/submm SEDs observed in large statistical surveys 
in a fast and flexible way, deriving physical parameters on an object-by-object basis. We also identify a non-parametric test of the 
fidelity of the model in practical applications through comparison of the model predictions for FIR colour and surface brightness 
with the corresponding observed quantities. This should be effective in identifying objects such as AGNs or star-forming galaxies 
with markedly different geometries to those of the calibrators of Xilouris et al. The results of the calculations are made available in 
the form of a large library of simulated dust emission SEDs spanning the whole parameter space of our model, together with the 
corresponding library of dust attenuation calculated using the same model. 

Key words. Radiative transfer - Scattering - (ISM): dust, extinction - Galaxies: spiral - Galaxies: individual: NGC 891 - Infrared: 
galaxies - submillimeter - Ultraviolet: galaxies 



1. Introduction 



Copious quantities of interstellar dust are present in the disks 

* We dedicate this paper to the memory of AngelosMisiriotis, sorely of all metal rich star-forming galaxies. This dust pervades all 
missed as a friend, collaborator and exceptional scientist. components of the interstellar medium (ISM), ranging from the 
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diffuse ionised and neutral medium filling most of the volume 
of the gaseous disk, through embedded neutral and molecular 
clouds of intermediate sizes and densities, down to the dense 
cloud cores on sub-parsec scales which are the sites of forma- 
tion of stars. The ubiquity and high abundance (relative to the 
available pool of interstellar metals) of grains is a fundamental 
result of the solid state being the favoured repository for refrac- 
tory elements in all but the coronal component of the ISM, and 
affects the very nature of galaxies. In particular, dust plays a ma- 
jor role in determining the thermodynamic balance of the ISM, 
through photoelectric heating and inelastic interactions with gas 
species, influencing the propensity of galaxies to accrete, cool 
and condense gas into stars (see e.g. Popescu & Tuffs 2010). In 
view of the physical role played by dust particles in the process 
of star formation, it can be regarded as a minor perversity of na- 
ture that, by virtue of their strong interaction cross section with 
stellar photons, the very same particles strongly inhibit and dis- 
tort our view of the resulting stellar populations, thus preventing 
a straightforward confrontation of theories for star formation in 
galaxies with observations. 

Statistical studies of the luminosity and colour distribu- 
tions of large samples of galaxies seen at different orientations 
(most notably the attenuation-inclination relation) have shown 
the severity of this effect even in the optical bands, and not just 
in exceptionally opaque systems, but in the major part of the 
population of star-forming disk galaxies in the local Universe 
(Driver et al. 2007, Shao et al. 2007, Choi 2007, Driver et al. 
2008, Unterborn & Ryden 2008, Padilla & Strauss 2008, Cho 
& Park 2009, Mailer et al. 2009, Ganda et al. 2009, Masters et 
al. 2010). This means that we not only need to solve the prob- 
lem of deriving intrinsic stellar emissivities from newly formed 
stars emitting predominantly in the ultraviolet, but also need to 
decode the effect of dust on the amplitude and colours of opti- 
cal light to derive the SF history of galaxies. This is fundamen- 
tally an ill-constrained problem using broad-band optical data 
alone, due to the age/metallicity/opacity degeneracy (Silva & 
Elston 1994, Worthey 1994, Li & Han 2008). Even for galax- 
ies bright enough for constraints on dust attenuation to be made 
through optical spectroscopy, for example from measurements 
of the emission from hydrogen recombination lines, large un- 
certainties still prevail due to the complex geometry of dust in 
galaxies and the different effect of this dust on the light from 
stellar populations of different ages. 

In recent years the advent of spaceborne infrared astronomy 
has meant that we have now the capability of viewing the ab- 
sorbed energy of the starlight in galaxies, which would seem to 
offer a solution to unravelling the effects of dust, at least from 
an observational perspective. Star-forming disk galaxies are ob- 
served to have most of this absorbed energy re-radiated in the 
far-infrared (FIR), with a significant amount also in the mid- 
infrared (MIR) and in the Polycyclic Aromatic Carbon (PAH) 
bands (see review of Sauvage et al. 2005). Simple energy bal- 
ance measurements indicate that indeed a large fraction of the 
stellar light is absorbed even in relatively quiescent spiral galax- 
ies in the local universe (Popescu & Tuffs 2002, Xu et al. 2006, 
Driver et al. 2008). Commonly, the total re-radiated energy is 
estimated by fitting the (typically sparsely sampled) measure- 
ments of PAH/dust emission either with phenomenological SED 
models (Devriendt et al. 1999, Dale & Helou 2002, Draine & Li 
2007) or with SED models that use templates (Xu et al. 1998, 
Sajina et al. 2006, Marshall et al. 2007, da Cunha et al. 2008). 
However, this approach is insufficient to derive the intrinsic UV 
and optical emission of galaxies, even in cases where the total 
infrared emission can be accurately estimated. 



There are several reasons why obtaining information about 
the relative absorption probabilities of UV/optical photons as a 
function of wavelength is a very complex problem, which can- 
not be solved by means of the application of single wavelength 
dependent correction factors as is the case for extinction correc- 
tions for stars Q. 

Firstly, as already noted, the ISM is highly structured, with 
dust present in both large scale and small scale structures. This 
strongly affects the dependence of attenuation on UV/optical 
wavelength, even for uniform geometrical distributions of stars, 
as shown by the study of attenuation of starlight in a turbu- 
lent ISM with a log-normal density distribution by Fischera & 
Dopita (2005). Similarly, various studies have shown that dif- 
ferent mass fractions of dust in clumps and diffuse structures 
can affect the amplitude and wavelength dependence both of the 
UV/optical attenuation and of the dust emission (Kuchinski et 
al. 1998, Bianchi et al. 2000a, Witt & Gordon 2000, Misselt et 
al. 2001, Misiriotis & Bianchi 2002, Pierini et al. 2004, Tuffs et 
al. 2004, Bianchi 2008). All these studies indicate that models 
for the attenuation of starlight and its re-emission in the infrared 
will only have predictive power if independent constraints on the 
geometry of the dust-bearing structures in the ISM are available. 

Secondly, stellar populations of different ages have system- 
atically different geometrical relations to these dust structures 
(Silva et al. 1998, Chariot & Fall 2000, Popescu et al. 2000) 
due to the disruption of birth clouds through the action of stel- 
lar winds and supernovae, the migration of stars away from these 
clouds with time and the subsequent dynamical evolution of stel- 
lar populations on Gyr timescales from a thin disk to a thick disk 
configuration (Wielen 1977). 

Thirdly, the response of grains to light is not just a question 
of the wavelength and flux of the photons hitting them but also of 
their optical properties (dependent on shape, size and composi- 
tion; see review by Draine 2010). A particularly complex aspect 
of this is the account that should in principle be made of any sys- 
tematic differences in composition on location in the ISM, such 
as for example the transition between pure refractory species in 
the diffuse ISM irradiated by UV to ices and ice-coated grains 
in opaque molecular clouds. Furthermore, even supposing the 
composition and size distribution is known, the stochastic emis- 
sion of small dust grains responding to impulsive heating has to 
be calculated throughout the volume, which is computationally 
challenging. 

Arguably, the most fundamental and challenging of these 
problems is the determination of the geometry of the dust- 
bearing structures in the ISM, since knowledge of this is a pre- 
requisite to obtaining a self-consistent solution embracing the 
amplitude and geometry of the stellar populations and the optical 
properties of the grains. Where FIR data is available, a powerful 
way of constraining the geometry of dust on all scales is to use 
the fact that grains act as test particles with FIR colours char- 
acteristic of the intensity and colour of the interstellar radiation 
fields (ISRF), which in turn are a strong function of morphologi- 
cal structures within galaxies. For example grains locally heated 
by the strong radiation fields within embedded star formation re- 
gions emit an infrared radiation that peaks around 60 /mi, while 
grains heated by the diffuse radiation fields in the disk of galax- 

1 We note that attenuation is an integral property of an extended dis- 
tribution (e.g. a galaxy) of light and should not be confused with the 
extinction measured for point sources (e.g. single stars). While the ex- 
tinction is simply proportional to the column density of dust and its 
wavelength dependence is determined by the optical properties of the 
grains, the attenuation depends on the distribution of dust and stars, the 
variation of dust properties and the orientation of the galaxies 



Cristina C. Popescu et al.: Modelling the spectral energy distribution of galaxies. 



3 



ies will emit an infrared radiation that typically peaks well long- 
wards of 100 /mi, but also exhibit substantial emission in the 
MIR region due to the stochastic heating of small grains (see 
Sauvage et al. 2005). In general using dust grains as tracers of 
the ISRF is a key way to constrain both opaque and translu- 
cent components of the ISM on all spatial scales, and the tech- 
nical implementation of this approach is to use radiative trans- 
fer (RT) calculations. Much effort has been put into this prob- 
lem of developing efficient RT codes, and there are now sev- 
eral available in the literature (see Kylafis & Xilouris 2005 for a 
comprehensive review on radiative transfer techniques used in 
galaxy modelling). Specific RT calculations of the integrated 
MIR/FIR/submm emission from grains heated by starlight in 
disk galaxies have been made for various specifications of the 
geometry of stars and dust by Siebenmorgen & Kriigel 1992, 
Silva et al. 1998, Popescu et al. 2000, Efstathiou & Rowan- 
Robinson 2003, Piovan et al. 2006, Bianchi 2008. 

Further empirical constraints on the geometry of dust in 
galaxies are possible for systems where major dust-bearing com- 
ponents of the ISM are translucent, thus allowing information 
about the distributions of dust and stars on scales of hundreds of 
parsecs to be inferred from images of the optical emission. This 
approach, which also relies on modelling of the data with RT 
techniques, was adopted by Xilouris et al. (1999) in their inves- 
tigation of the large scale diffuse components in edge-on spiral 
galaxies, revealing reproducible trends in the geometry of stars 
with respect to the dust. The analysis of Xilouris et al. was also 
notable in that it was also able to measure the extinction law for 
the visible component of dust for these systems, showing it to 
be consistent with that of dust in the diffuse ISM of the Milky 
Way in the observed optical/NIR range. In principle the geomet- 
ric constraints on the distributions of stellar emissivity and dust 
derived from optical observations of translucent regions can be 
extended to more obscured components of the ISM, such as for 
example the molecular layer in the disk out of which stars are 
forming, by making use of the physical connection between gas, 
stars and dust to further constrain the geometry of the problem. 
Overall, our knowledge of galaxies, at least in the local Universe, 
suggests that the specification of geometry does not need to be 
completely ad hoc and can be empirically constrained by the 
colour of the integrated dust emission, by the optical measure- 
ments of structure, and by considerations of physical plausibility. 
This approach can be directly applied to the panchromatic obser- 
vations to decode their information into basic physical parame- 
ters of galaxies, and we therefore call this approach decoding 
observed panchromatic SEDs (see Popescu & Tuffs 2010). 

The only alternative to invoking empirical constraints on the 
distributions of dust and stars in galaxies is to calculate these 
geometry from first principles, using numerical simulations of 
how galaxies form and evolve. This approach was followed 
by Chakrabarti et al. (2008), Chakrabarti & Whitney (2009), 
Jonsson et al. (2009) and we call it encoding predicted physi- 
cal quantities. A particular advantage of the encoding approach 
is in its application to high redshift galaxies with more complex 
geometries due to frequent mergers and interactions, where hy- 
drodynamical simulations are ideal for describing structures on 
scales greater than ca. 100 pc (Jonsson et al. 2009), at least in a 
statistical sense. The coupling of cosmological simulations with 
SED modelling tools is moreover an obvious way to interpret the 
cosmological evolution of luminosity functions of both the direct 
and re-radiated components of stellar light, something which has 
already been done using semi-analytic models of structure for- 
mation (Almeida et al. 2010). On the other hand, the encoding 
process relies on the assumption that the theory used to predict 



the relation between gas, dust and stars in galaxies is complete. 
Furthermore, deriving information on galaxies on an object-by- 
object base is better suited to the decoding process, at least for 
evolved disk galaxies in the local Universe with more uniform 
geometries. 

In this paper we follow the decoding approach, using an 
empirical determination of geometry to self-consistently pre- 
dict the UV/optical attenuation of starlight and the correspond- 
ing dust/PAH emission SEDs of star-forming galaxies in the 
local Universe. To this end we use an updated and enhanced 
version of the original model of Popescu et al. (2000, here- 
after Paper I) which incorporates the geometrical constraints on 
dust and stars in spiral galaxies found by Xilouris et al. (1999). 
Specifically, we calculate a comprehensive library of spatially 
integrated dust/PAH re-emission SEDsEI of disk galaxies as a 
function of a minimal set of physical parameters needed to pre- 
dict the dust/PAH emission. We further describe how this set of 
dust/PAH re-emission SEDs can be self-consistently combined 
with a new library of UV/optical dust attenuations calculated 
using the same model (an update of the existing library from 
Tuffs et al. 2004; hereafter Paper III) to invert an observed set 
of broad-band photometry of a galaxy spanning the UV/optical 
- FIR/submm range to derive the intrinsic (i.e. as would be ob- 
served in the absence of dust) UV/optical emission of the galaxy. 
In particular we demonstrate that this analysis can, without sig- 
nificant loss in accuracy, be done independent of a priori as- 
sumptions of the detailed mathematical form of the dependence 
of SFR on time. This approach allows the UV/optical SED to be 
dereddened in a fast and flexible way without recourse to popula- 
tion synthesis models. The latter models can then be compared to 
the dereddened UV/optical SED to derive the SF history in a fur- 
ther, independent step, avoiding bias due to the age/metallicity- 
opacity degeneracy. 

The primary motivation to develop this approach was to ex- 
tend the applicability of a fully self-consistent RT solution to 
the large statistical samples of optically selected local Universe 
star-forming galaxies spanning a full range of mass, morphology 
and environment for which, thanks to facilities such as GALEX, 
WISE, Spitzer, AKARI and Herschel (e.g. Driver et al. 2009, 
Eales et al. 2010, Martin 2010) integrated photometry is now for 
the first time becoming routinely available across the full UV- 
submm ranged In addition to showing how the UV/optical atten- 
uation is constrained through the observed colour and amplitude 
of the observed PAH/dust emission spectrum, we also describe 
how to utilise morphological information from high resolution 
optical observations of galaxies (such as linear sizes of disks 
and the bulge-disk decompositions) which are expected to come 
from the next generation of wide field imaging spectroscopic 
surveys of local Universe galaxies (e.g. Driver et al. 2009). In 



2 All available in electronic format at CDS data base 

3 Here we emphasise once more that our model has been calibrated 
for local Universe galaxies, and therefore its applicability should mainly 
lie within low redshift galaxies. The models are also targeted to spiral 
galaxies. We have not attempted to model elliptical galaxies, starburst 
galaxies or AGN nuclei. Models for starburst galaxies can be found in 
Rowan-Robinson & Efstathiou 1993, Kriigel & Siebenmorgen 1994, 
Silva et al. 1998, Efstathiou et al. 2000, Takagi et al. 2003, Dopita et 
al. 2005, Siebenmorgen & Kriigel 2007, Groves et al. 2008. Models 
for starburst dwarfs can be found in Galliano et al. (2003). Models 
for AGN torus were presented by Pier & Krolik 1992, Efstathiou & 
Rowan-Robinson 1990, Granato & Danese 1994, Efstathiou & Rowan- 
Robinson 1995, Nenkova et al. 2002, Dullemond & van Bemmel 2005, 
Fritz et al. 2006, Honig et al. 2006, Schartmann et al. 2008, Nenkova et 
al. 2008. 
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keeping with the decoding approach that we have adopted, we 
envisage that the model will be well suited to the establishment 
of empirical relations between SF activity and SF history of op- 
tically selected galaxies with a large range of properties of in- 
terest, such as dynamical mass (both of the parent DM halo and 
in disk), baryonic gas content, dust content, relativistic particle 
content, and environment. Amongst other applications the model 
should also be useful for evaluating the relation between broad- 
band and spectroscopic measures of both SFR and opacity in 
spiral galaxies so that these quantities may be measured in cases 
where a full panchromatic coverage of the UV-submm SED is 
not yet available. 

Although we have calibrated the model using the optical ob- 
servations of the galaxies from Xilouris et al. (1999) which are 
all spirals in the middle range of the Hubble sequence, we have 
incorporated a bulge-to-disk luminosity ratio parameter which 
will extend the applicability of the model to galaxies of earlier 
and later morphological types. We have also identified a non- 
parametric test to evaluate the fidelity of this model by virtue 
of its prediction of the colour-surface-brightness relation for spi- 
ral galaxies. This test will allow galaxies with systematic differ- 
ences from the geometry adopted in this model to be identified. 

This paper is organised as follows. In Sect. [2] we describe 
our model for the dust emission, together with the physical mo- 
tivation associated with its underlying assumptions, and for- 
mulate the mathematical framework underpinning the quanti- 
tative comparison of dust emission with the attenuation calcu- 
lations presented in this paper. We draw particular attention to 
the changes in the model compared to the original formula- 
tion in Paper I. Specifically, we describe the inclusion of PAH 
molecules, needed to extend the applicability of the model to 
shorter (mid-IR) wavelengths, and which also led to a revision 
to the overall grain size distribution to preserve consistency with 
the extinction and emission properties of interstellar dust. A fur- 
ther major update to the model is the inclusion of the contri- 
bution to dust heating of the young stellar population at optical 
wavelengths, avoiding the artificial cut-off of the emission of this 
population beyond the UV range in the original version of the 
model from Paper I. In this section we further describe a revised 
treatment of the emission from star-formation regions, incor- 
porating a wavelength-dependence of the escaping UV light to 
achieve consistency with their attenuation properties as formu- 
lated in Paper III and including an updated template for the cor- 
responding infrared emission using the recent model of Groves 
et al. (2008). In Sect. [3] we illustrate our model for the case of 
the well studied galaxy NGC 891. A more general quantitative 
analysis of the effects of the main assumptions of the model on 
the predicted spatially integrated SEDs as well as tests on the 
fidelity of the model are given in Sect. |4] We then present our 
library of simulated dust emission SEDs in Sect. [5] and discuss 
the predicted variation of these SEDs with the main parameters 
of the model in Sect. [6] We summarize the main results of the 
paper in Sect. [8] 

The updated calculations of attenuation of stellar light for 
the diffuse component are given in AppendixlBl To facilitate the 
use of this library of attenuation in combination with the dust 
emission calculations we give in Appendix [C] formulae for the 
composite attenuation of stellar light from the different geomet- 
rical components of stellar emissivity in terms of the parameters 
used for the prediction of the dust emission. For an easy access 
to potential users of the model we give in Appendix|D]a step-by- 
step guide on how to use the model for fitting real data, together 
with the full mathematical formalism. In this appendix we also 
describe the non-parametric test to objectively judge the efficacy 



of the model by virtue of its prediction of the colour-luminosity 
relation for spiral galaxies. 

2. Calculation of the infrared SEDs 

We calculate the infrared SEDs using an updated version of the 
model from Paper I which was applied to the interpretation of 
dust emission from edge-on galaxies by Misiriotis et al. (2001; 
Paper II) and also used to make predictions for the effects of 
dust in the optical range in Paper III and by Mollenhoff et al. 
(2006; Paper IV). In the following we describe the changes in 
detail, giving their physical motivation, and showing how they 
interlink with the mathematical formulation of the model. 

2. 1 . The distribution of stars and dust 

Star-forming galaxies are fundamentally inhomogeneous, con- 
taining highly obscured massive star-formation regions, as well 
as more extended large scale distributions of stars and dust. It 
is important to emphasize that the large-scale distribution of dif- 
fuse dust plays a major role in mediating the propagation of pho- 
tons in galaxy disks and dominates the total bolometric output of 
dust emission. The discovery of this diffuse component was one 
of the highlights of the Infrared Space Observatory (ISO)(see 
Tuffs & Popescu 2006 and Sauvage et al. 2006 for reviews on 
the ISO science legacy on normal nearby galaxies). This result 
was obtained both from analysing the integrated properties of in- 
frared emission from galaxies (Popescu et al. 2002) but also from 
resolved studies of nearby galaxies (Haas et al. 1998, Hippelein 
et al. 2003, Tuffs & Gabriel 2003, Popescu et al. 2005) and is 
now being confirmed by the infrared data from Spitzer (Perez- 
Gonzalez et al. 2006, Hinz et al. 2006, Dale et al. 2007, Bendo 
et al. 2008, Kennicutt et al. 2009) and AKARI (Suzuki et al. 
2007). Our model being empirically motivated we tried to in- 
corporate all available observational constrains provided by the 
data. Accordingly our model includes both a diffuse component 
and a localised component representing the star-formation com- 
plexes. 

A schematic picture of the geometrical components of the 
model is given in Fig.[T] together with a mathematical prescrip- 
tion of the stellar emissivities and dust opacities used in the 
model. 

The large scale distribution of stars and dust are approxi- 
mated as continuous spatial functions of stellar emissivity and 
dust opacity, which we refer to as "diffuse" distributions. We 
have separate distributions for the old and young stellar popula- 
tions, and we also consider separate distributions for diffuse dust 
associated with these populations. 

The old stellar population resides in a disk and a bulge, 
with its emissivity described by a double exponential and a de 
Vaucouleurs distribution, respectively: 

^,^ Z )^^,0,0)exp(-^-|L) 

+ 77 bulge (^, 0, 0) exp(-7.67 B 1/4 ) fi~ 7/8 , (1) 
= ^R 2 + z 2 (a/b) 2 

Re 

where R and z are the cylindrical coordinates, 77 dlsk (/l, 0, 0) is the 
stellar emissivity at the centre of the disk, hf sk , zf sk are the scale- 
length and scaleheight of the disk, rj hulge (A, 0, 0) is the stellar 
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Fig. 1. Schematic representation of the geometrical distributions of stellar and dust emissivity together with a mathematical prescrip- 
tion of the stellar emissivities and dust opacities used in the model. Here, and in the main body of the text we use the superscript 
"disk", "bulge" and "tdisk" for all the quantities respectively describing the disk (the old stellar disk plus the associated dust disk, 
whether the latter is sometimes referred to as the "first dust disk"), the bulge and the thin disk (the young stellar disk plus the 
associated dust disk, whether the latter is sometimes referred to as the "second dust disk"). 



emissivity at the centre of the bulge, R e is the effective radius 
of the bulge, and a and b are the semi-major and semi-minor 
axes of the bulge. 

The total luminosity of the old stellar population in the disk 
is then given by: 



r disk 



:4^77 disk (i,0,0)zf Sk (/!f Sk ) 2 



(3) 



whereas for the bulge there is no exact analytical formula for the 
spatially integrated luminosity. 

The dust associated with the old stellar population ("the first 
dust disk") is also described by an exponential disk: 



k?£(A,R,z) = CO*. 0,0) exp 



R 

/disk 
"d 



-disk 



(4) 



where K^(A, 0, 0) is the extinction coefficient at the centre of 
the disk and h dlsk and z* sk are the scalelength and scaleheight of 
the dust associated with the old stellar disk. 

In our model the parameters describing the geometry of the 
old stellar population and of the first dust disk are empirically 
constrained from resolved optical and near-IR images via the re- 
sults of the modelling procedure of Xilouris et al. (1999). For 
edge-on systems these calculations completely determine the 
scale heights and lengths of the exponential disk of old stars and 
associated diffuse dust, as well as the effective radius and ellip- 
ticity of the dustless stellar bulge. This is feasible for edge-on 
systems since the scale height of the dust is less than that of the 
stars. By analysing a sample of 5 nearby edge-on galaxies with 
morphological classification in the range of Sb-Sc, Xilouris et 
al. found that the scalelength of the dust was larger than that of 
the stellar disk, that the scaleheight of the stars was larger than 
that of the dust and that the scalelength of the stellar disk de- 
creases with increasing wavelength. In this way Xilouris et al. 



was also able to find a general relation between the scaleheights 
and scalelength of old stars and dust and the dependence of this 
relation on wavelength. The derived relation enabled us to fix 
the relative distribution of old stars and dust in our model. In 
Papers I and II we verified that the derived geometrical param- 
eters from Xilouris et al. correctly predicted the dust emission 
SEDs of the 5 nearby edge-on galaxies. Observationally it was 
also confirmed that the dust disk not only has a larger scalelength 
than the stellar disk (Alton et al. 1998, Davies et al. 1999), but 
was detected to physically extend well beyond the stellar disk 
(Popescu & Tuffs 2003; see also Popescu et al. 2002 and Hinz 
et al. 2006 for dwarf galaxies). From resolved studies of galax- 
ies it was also found that there is a large scale distribution of 
diffuse dust having a face-on opacity that decreases with radius 
(Boissier et al. 2004, Popescu et al. 2005, Perez-Gonzalez et al. 
2006, Boissier et al. 2007, Munoz-Mateos et al. 2009). 

The geometrical parameters of our model are listed in 
Table IE. ll where all the length parameters are normalised to the 
B-band scalelength of the disk, hf sk (B). In our calculations we 
take hf sk (B) = hj^j = 5670 pc, the fixed reference scalelength 
of our model galaxy, as derived for NGC 891. 

The "young" stellar population and associated dust are also 
specified by exponential disks, which are taken to have small 
scaleheights (thin disks, thereby the superscript 'tdisk'), which 
we shall refer to as the "young stellar disk" and the "second dust 
disk": 



T] mk (A, R,z) = r] ldisk (A, 0, 0) exp I - 



R 



.tdisk 



K^ k (A,R,z) = K^ sk 0l,O,0) exp 



R 

/tdisk 
"d 



\z\ 

-tdisk 



\z\ 

.tdisk 



(5) 



(6) 
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where if dmk (A, 0, 0) is the stellar emissivity at the centre of the 
thin disk, hf lsk and z tdlsk are the scalelength and scaleheight of 
the thin disk, k^\A, 0, 0) is the extinction coefficient at the cen- 
tre of the thin disk and h tdlsk 



and zf sk 



are the scalelength and 
scaleheight of the dust associated with the young stellar disk. 
The total luminosity of the young stellar disk is given by: 



Lf sk = 4 7T ij ms \A, 0, 0) zf sk (hf sk ) 2 



(7) 



Unlike the old components, the parameters describing the 
geometry of the young components cannot be constrained from 
images of stellar light, because the young stellar populations are 
highly obscured in most cases. It was thus necessary to constrain 
these parameters from physical considerations. The scale height 
of the young stars was taken to be 90 pc (the value for the Milky 
Way) and its scale length is equated to that of the "old stellar 
disk" in B-band, /z tdisk = h disk (B). The dust associated with the 
young stellar population was fixed to have the same scalelength 
and scaleheight as for the young stellar disk. The reason for this 
choice is that our thin disk of dust was introduced to mimic the 
diffuse component of dust which pervades the spiral arms, and 
which occupies approximately the same volume as that occupied 
by the young stars. This choice is also physically plausible, since 
the star-formation rate is closely connected to the gas surface 
density in the spiral arms, and this gas bears the grains which 
caused the obscuration. In principle it would be more realistic 
to place the young stellar population and associated dust into a 
spiral pattern rather than into a disk. However, in practice, as we 
show in Sect. 14. 21 this makes little difference to the predicted vol- 
ume integrated dust emission SED and starlight attenuation of a 
galaxy. Furthermore, the inclusion of spiral arms would neces- 
sitate another model parameter to describe this more complex 
geometry. Also, in most practical applications for the interpre- 
tation of the integrated emission of galaxies in large statistical 
samples, photometric information about the spiral arm pattern is 
not available. For all these reasons we consider the approxima- 
tion of the second dust disk to be both reliable and practical. The 
geometrical parameters of the young stellar disk and second-dust 
disk are listed in Table IE. ll 

Apart from the geometrical parameters described above, the 
distributions of diffuse stellar emissivity and dust are also de- 
scribed in terms of their amplitudes. The amplitudes of the two 



dust disks , can be expressed in terms of the central 
face-on opacity in the B band, T ^ d,sk ,t^' lsk , defined by: 



t^*** =2 K f£(A B ,0,0)zf sk 



f, tdisk _ ~ tdisk/ i n n\ tdisk 
R = 2 «ext 0**»0,U)Z d 



(8) 



(9) 



The opacity of the first dust disk can be derived jointly with 
the geometrical parameters from the optimisation technique of 
Xilouris et al. (1999). The optimisation also determines the ex- 
tinction law of the diffuse dust empirically, since the calcula- 
tions are done independently for each optical/NIR image. For 
the galaxies studied by Xilouris the derived extinction law was 
consistent with a Milky Way type dust, which is also the type 
of dust adopted for our dust model (see Sect. 12.41 and lAb. The 
opacity of the second dust disk was a free parameter in the cal- 
culations from Papers I and II, and is strongly constrained by the 
level of the submm emission. To minimise the number of free 
parameters, we fix 



.disk 



f, tdisk 



(10) 



to the value 0.387 found for our proto-type galaxy NGC 891, 
which is also close to what was found for a second edge-on 
galaxy with submm data, NGC 5907, which we modeled in 
Paper II. We note here that the attenuation-inclination relation 
predicted for this fixed ratio of opacities in the two dust disks 
was found to successfully reproduce the observed attenuation- 
inclination relation of a large and statistically complete sam- 
ple of galaxies from the Millennium Galaxy Catalogue Survey 
(Driver et al. 2007). We thus adopt as a free parameter of the 
model the total central face-on opacity in the B-band r 



J. 



f f.disk fjdisk 

t b = r D + ti 



(11) 



The remaining parameters - the amplitudes of the two stellar 
disks 77 disk , 77 tdisk and of the bulge rj bul s e - and their link to the free 
parameters of the model are discussed in Sect. 12.31 

2.2. The clumpy component 

An important component of the spatially integrated dust emis- 
sion from star-forming galaxies arises from dust in the birth- 
clouds of massive stars, as previously modelled by Silva et al. 
(1998), Popescu et al. (2000), Chariot & Fall (2000), Efstathiou 
& Rowan-Robinson (2003), Jonsson et al. (2009). Because these 
clouds are spatially correlated with their progeny on parsec 
scales, they are illuminated by a strong UV-dominated radiation 
field of intensity 10-100 times that in the diffuse ISM. This 
gives rise to a localised component of emission from grains in 
thermal equilibrium with these intense radiation fields, which, 
despite the tiny filling factor of the SF regions in the galaxy 
can nevertheless exceed the entire diffuse infrared emission of 
a galaxy at intermediate wavelengths (ca. 20 to ca. 60 //m). It is 
therefore particularly important to incorporate a clumpy compo- 
nent of dust associated with the opaque parent molecular clouds 
of massive stars. Following Popescu & Tuffs (2005) we refer 
to these clumps as "active clumps". The active clumps are as- 
sumed to have the same spatial distribution as the young stellar 
disk and the second dust disk. Furthermore, it is assumed that 
the properties of these clumps do not systematically depend on 
their radial location within the galaxy. In reality we expect star 
formation complexes in more pressurised regions (such as the 
inner disk) to be more compact and therefore have warmer FIR 
colours than their counterparts in low pressure regions (such as 
the outer disk), as modelled by Dopita et al. (2005). And in- 
deed this phenomenon has been observed in HII regions in M 33 
(Hippelein et al. 2003). However, bearing in mind that we will 
empirically calibrate our template (see Sect. 12.81 1 on the whole 
ensemble of star forming complexes in the Milky Way, the as- 
sumption we make in this regard should not significantly affect 
our predictions for the integrated emission of galaxies. 

Since birthclouds of stars are typically fragmented due to the 
combined effects of supernovae, stellar winds, and the general 
motion of stars away from the clouds, only a certain fraction of 
the total luminosity of massive stars in a galaxy will be locally 
absorbed. Following the original formulation in Paper I we de- 
note this fraction in our model by the "dumpiness factor", F. 
Since birthclouds are completely opaque at all UV/optical/near- 
IR wavelengths (e.g. Sievers et al. 1991), F can be physically 
identified with the luminosity-weighted mean fraction of direc- 
tions from the massive stars, averaged over the lifetime of the 
stars, which intersect the birthcloud. This concept, first intro- 
duced by Popescu et al. (2000), allows UV light to freely stream 
away from star-forming regions in some fraction of directions, 
allowing the diffuse dust to be illuminated with more UV pho- 
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tons than would have been the case if the young stars had been 
assumed to be completely cocooned in their parent molecular 
clouds. This is in qualitative accordance both with the high fre- 
quency at which counterparts of star-formation regions seen at 
24 fim are detected in the UV, as well as with the predominance 
of PAH emission from diffuse dust in galaxies correlated with 
cold dust (Bendo et al. 2008). 

A more detailed description of the clumpy component is 
given in Sect. l2.5.Tl and l2.8l 

An approximation of our model is that the heating of the 
grains in the active clumps is dominated by photons from the 
stellar progeny, and we neglect any external contribution from 
the diffuse ambient radiation fields in the galaxy. This should 
be an excellent approximation for spiral galaxies, where the fill- 
ing factor of star formation regions is small. We note that if this 
approximation were invalid, we would be forced to perform ra- 
diative transfer calculations with parsec resolution to properly 
describe the heating of grains in the optically thick birth clouds, 
which would render the calculation of a comprehensive set of 
SEDs as given in this paper intractable with current computing 
resources. To date, the best resolution achieved with an adap- 
tive grid code when modeling the observed dust emission SED 
from an individual galaxy is about 20 pc (Bianchi 2008). In any 
case, we would only expect collective effects, whereby photons 
from adjacent SF regions provide a significant fraction of the 
dust heating, to be significant in galaxies with very high volume 
densities of SF regions, such as the central region of a starburst 
galaxy. The low filling factors of opaque clouds in spiral disks is 
supported by high resolution surveys of the Galactic Plane with 
large ground-based telescopes in the submm (e.g. the APEX 
Telescope large Area Survey - ATLAS; Schuller et al. 2009). 
At these wavelengths the clouds are optically thin, thus directly 
tracing the total column density of dust, yet only sporadic peaks 
of dust emission can be seen with local enhancement in visual 
optical depth much greater than 1 . This contrasts with surveys 
in the CO lines (Matsunaga et al. 2001), which show a more 
highly filled distribution of emission due to optical depth effects 
in the radio molecular line. Recently the new Herschel Infrared 
Galactic Plane Survey (Hi-Gal) (Molinari et al. 2010) has re- 
vealed the fainter, more extended emission components which 
dominate the morphology of the dust emission in the Milky Way, 
showing again that the optically thick molecular cores have only 
a small filling factor. In summary, based on all available obser- 
vational evidence we believe it is a reasonable approximation to 
ignore the external heating of active clumps. 

Another approximation of our model is that the clumpy dis- 
tribution of dust is exclusively associated with the opaque parent 
molecular clouds of massive stars - the active clumps. In reality, 
some of the dust in the diffuse dust disks may also be in clumps 
without internal photon sources, and thus heated only by the dif- 
fuse ambient radiation fields in the disk. Following Paper I (see 
also Popescu & Tuffs 2005) we refer to such clumps as "pas- 
sive" or "quiescent" clumps. Provided passive clumps are op- 
tically thin, the transfer of radiation through the disks will be 
virtually identical to that in a homogeneous disk. However, once 
the passive clumps become optically thick, the self-shielding of 
grains will yield a solution with a reduction in both the attenu- 
ation of the stellar light and the infrared emission compared to 
the case where the same mass of grains is diffusely distributed. 
These effects have been quantified by Bianchi et al. (2000a) who 
showed that the shape of the infrared SED of a disk galaxy can 
be strongly affected shortwards of 200 fim (see also Misselt et 
al. 2001 for passive clumps distributed in a spherical shell ge- 
ometry). We do not include passive clumps in our model for 



both physical and empirical reasons. Physically, it seems likely 
that most passive clumps will be optically thin, since as soon as 
they become optically thick to the impinging external UV light, 
they will lose their principal source of heating (the photoelec- 
tric effect) and will be prone to collapse and form stars (Fischera 
& Dopita 2008). Empirically, this presumption is supported by 
the recent findings of Holwerda et al. (2007a,b), that the struc- 
ture of the diffuse ISM consists of optically thin dusty clouds. 
Furthermore, the change in the shape of the predicted dust emis- 
sion SEDs imposed by the incorporation of passive clumps tends 
to provide a colder solution than needed to fit real data (Bianchi 
2008). 

2.3. Constraints on the intrinsic SEDs of the stellar 
populations 

As mentioned in Sect. 12.11 the amplitudes of the stellar popu- 
lations of the two stellar disks (Tj Aisk (A, 0, 0) and ?7 tdisk (/i, 0, 0)) 
are parameters that are yet to be determined. In other words we 
want to find a solution for the attenuation of stellar light and dust 
emission SED that can fit observed SEDs to provide resulting 
intrinsic stellar SEDs. At the same time, in our modelling pro- 
cedure we deliberately do not want to use population synthesis 
models to fully fix the intrinsic SED, as we want to be as free as 
possible of any assumptions about the SF history. Nonetheless, 
some constraints on the input SEDs need to be made to avoid 
having the amplitude of the stellar SED a free parameter at each 
of the sampled wavelengths listed in Table IE. 21 Since the main 
factors shaping the dust emission SEDs are the total luminosity 
of the young stellar population and of the old stellar populations, 
we choose to have these two quantities as free parameters of our 
model and produce calculations for all combinations of these 
two variables. We then assume that each of the two components 
(young and old stellar populations) have a fixed wavelength de- 
pendence (a fixed template SED), thus reducing the number of 
free parameters describing the intrinsic SEDs of the stellar pop- 
ulations to two. We will show in Sect. l4.3l that these assumptions 
have a negligible effect on the predictions for the dust emission 
SEDs and the associated parameters t b and F. In turn, the solu- 
tion for the wavelength dependence of the attenuation in the UV- 
optical range, which depends on r^, F (and inclination), will be 
similarly secure, thus allowing the true dereddened stellar emis- 
sion SED to be recovered. In effect we are using the fact that 
the dependence of the dust emission SED on the colour of the 
stellar photon field depends primarily on the ratio between the 
luminosities of the young and old stellar populations rather than 
on the detailed colour of the emissions from either of these pop- 
ulations. This way of constructing the model will allow SF his- 
tories to be extracted from the dereddened stellar emission SED 
in a separate step from the extraction of r B and F. 

In constructing the template SEDs of the stellar populations 
we consider the term "optical", "UV", "ionising UV" and "non- 
ionising UV" to denote the wavelength ranges A > 4430 A, A < 
4430 A, A < 912 A, 912 < A < 4430 A, respectively, thereby 
marking the boundary between the UV and optical regime in the 
Bband. 

2.3.1 . The template SED of the old stellar population in the 
disk 

In defining the fixed shape of the SED of the old stellar pop- 
ulation in the disk we only consider optical radiation and ne- 
glect any contribution in the UV. The precise wavelength depen- 
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dence of this SED was fixed to the empirical relation obtained 
by Xilouris et al. (1999) from fitting simulated images produced 
by radiative transfer calculations to the observed B, V, I, J, K 
imaged of NGC 891. The SED of the old stellar population in 
the disk, as "calibrated" on NGC 891 is given as the light blue 
curve in Fig. [8] In applications to other galaxies this curve needs 
to be scaled according to the total output of the old stellar pop- 
ulation in the disk of the modelled galaxy. For this purpose we 
use a unitless parameter old defined as 



old — 



r disk 



-old 



(12) 



where L°~? = 2.241 x 10 37 W, which corresponds to 10 times 
the luminosity of the non-ionising UV photons produced by a 
lM Q /yr young stellar population, as defined in Sect. 12.3.31 Thus 
old is the normalised luminosity of the old stellar populations in 
the disk and is adopted to be one of the free parameters of the 
model. In this scheme the best fit solution for NGC 891 corre- 
sponds to old = 0.792. 

The intrinsic spectral luminosity densities of the old stellar 
population in the disk corresponding to old = 1, L° it , used in 



our calculations are listed in Table IE. 21 The spectral luminosity 
density of the old stellar population in the disk of a model galaxy 
is then given by: 



rdisk _ j j y j old 



(13) 



where L d v ' sk has been previously defined in terms of r\ Al ^(A, 0, 0) 
inEq.fJ 



,old 



opt 



;dA 

A, unit 



(14) 



and the value of L°~? can be derived by integrating over wave- 
lengtlflthe values of L°' d , from TablelR2l 

to LJ v, unit ' ' 

As already mentioned, we will show in Sect. 14.31 that, al- 
though the template SED of the old stellar population used in 
the model is fixed to a single shape (here defined empirically on 
the galaxy NGC 891), it is still valid to use it to model the dust 
emission SEDs of other galaxies which may have very different 
stellar emission SEDs. 



2.3.2. The template SED of the old stellar population in the 
bulge 

The wavelength dependence of the stellar luminosity produced 
by the bulge stellar population is simply linked to that of the 
old stellar population in the disk via a wavelength independent 
bulge-to-disk ratio B/D, one of the free parameters in our model. 



L b v ulge = (B/D)xLf sk 



(15) 



4 In the K and H bands we rescaled the luminosities derived by 
Xilouris et al. (1998) to more recent, higher quality data from the 
2MASS survey. 

3 Throughout the paper all spectral integrations are done over wave- 
length rather than over frequency. This ensures self-consistency in cases 
where functions are sparsely sampled in wavelength. 



2.3.3. The template SED of the young stellar population 

The wavelength dependence of the stellar luminosity produced 
by the young stellar population in the UV cannot be constrained 
empirically, as this population is heavily obscured by dust. 
Because of this we used population synthesis models. We used 
the models from Kotulla et al. (2009), making standard assump- 
tions for the star-formation history, namely an exponentially de- 
clining star-formation rate with a time constant r = 5 Gyr, so- 
lar metallicity and a Salpeter IMF with an upper mass cut-off of 
100 M . We emphasise that, although arbitrary, this choice of pa- 
rameters will not significantly bias the model dereddening of the 
observed UV/optical SEDs, since, as we will show in Sect. 14.31 
using different assumed shapes for the stellar emissivity SED 
does not significantly affect the predictions for the dust emission 
SEDs. 

We also consider that a fraction fi „- U v - 0.20 of the ionising 
UV photons (emitting shortwards of 912 A) from the massive 
stars within the HII regions are locally absorbed by dust. This 
makes only a minor contribution to the dust heating, an order 
of magnitude less than the contribution of the non-ionising UV 
photons. 

A new feature of the model is the consideration of the opti- 
cal emission from the young stellar population embedded in the 
second dust-disk. Previously the emissivity function of this pop- 
ulation had been artificially truncated longwards of the B band. 
In order to fix the shape of the SED template for the young stel- 
lar population in the optical we again make use of the best fit 
solution to NGC 891. Specifically, the optical radiation from 
this template at each wavelength sampling point was fixed to 
be the residual between the prediction of the population synthe- 
sis model for the best fit of NGC 891 and the empirical SED of 
the old stellar population (derived from the best fit of NGC 891 
as explained above j3 The resulting optical part of the SED tem- 
plate of the young stellar population was then combined with the 
UV part (which is simply the population synthesis SED for the 
best fit SFR for NGC 891) to fix the shape of the template over 
the full UV/optical range. For use in other model galaxies, the 
overall amplitude of this SED of the young stellar population, as 
"calibrated" on NGC 891, needs to be scaled according to the 
S FR of that galaxy, adopting S FR as a free parameter of the 
model. 

For logistical purposes only, we chose as a unit for the lumi- 
nosity of the young stellar population L?°"" 8 = 4.235 x 10 36 W, 
which is the luminosity produced by a SFR of 1 M Q /yr (for 
the standard assumptions about the SF history described be- 
fore) in the range (912. - 50000.) A. L; °™ s can be given as a 
sum of the luminosity of the non-ionising UV photons L^" 8 = 
2.241 x 10 36 W and the luminosity of the optical photons emitted 
by the young stellar disk L y ^ opt = 1-994 x 10 36 W. As already 
introduced in Sect. 12.3.11 we also defined the unit for the lu- 
minosity of the old stellar population L°~l to correspond to 10 



6 The fraction of ionising UV photons that is absorbed by dust in HII 
regions exhibits a broad range of values, varying from 0.3-0.7 (Inoue et 
al. 2001; Inoue 2001). However, even if this fraction approaches unity, 
their contribution to the dust emission of the star-forming complexes 
is still only at the percent level, because the intrinsic luminosity of the 
ionising UV photons is so much smaller than that of the non-ionising 
UV photons (Bruzual & Chariot 1993). This means their contribution 
to the total dust emission is even less. 

7 Physically, this approach corresponds to the dynamical heating of 
stellar populations born in the young stellar disk and their transfer over 
time into the larger scale height old stellar disk due to inelastic scatter- 
ings with spiral arms and/or giant molecular clouds. 
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times the luminosity of the non-ionising UV photons produced 
by a 1 Mo/yr young stellar population. In addition we also define 
L young _ Q 267 1Q 36 w f th i uminosit y f the ionising 

unit, wn-uv J to 

UV photons. We note that L^°"."f„ is not included in the defi- 



nition of L >ou " g . 

unit 



unit, wn—uv 



We then parameterised the calculations in terms of S FR 



SFR 
1 M yr 



r tdisk 



(16) 



where the SFR is linked to the normalised luminosity of the 
young stellar population in the thin disk. 

In an equivalent way we can link the normalised luminosity 
of the young stellar population in the UV to the SFR: 



SFR 

1 M yr" 



j mi 



tdisk 



, young 



(17) 



where L y °T g = 2.241 x 10 36 W, as defined previously, and 



SFR 



V. 



tdisk 



IMoyr- 1 L 



young 
unit, ion—uv 



(18) 



The intrinsic spectral luminosity densities of the young stel- 
lar population corresponding to S FR - 1 M yr _1 , L y °™ g , used 
in our calculations are listed in Table IE. 2l 

The spectral luminosity density of the young stellar popula- 
tion in the thin disk of a model galaxy is then given by: 



r tdisk 



■ young 



SFR r 

lM yr-' v ' unit 

r tdisk 



(19) 



where L|f has been previously defined in terms of T] (A, 0, 0) 
inEq.Hl 



jjoung 
unit 



I 



UV+opt 



A, unit 



(20) 



and the value of L y ° n "" 8 
length the values of L yom8 . from Table lE.2l 

fc v, unit ' 1 



2.4. The dust model 



can be derived by integrating over wave- 



The dust model - i.e. the prescription for the optical properties 
of the grains, the chemical composition and the grain size dis- 
tribution - was updated to include PAH molecules, which are 
regarded as the carriers of the unidentified infrared emission fea- 
tures commonly seen in the mid-IR emission spectrum of star- 
forming galaxies. As shown by Zubko, Dwek & Arendt (2004), 
there is no unique model incorporating PAHs that can simul- 
taneously fit the main observational constraints on such mod- 
els, namely the extinction and emission properties of the dif- 
fuse cirrus in the Milky Way. Several models have been pro- 
posed in the literature, including the models from Zubko et al. 
(2004), the model of Fischera & Dopita (2008) which consid- 
ers an extra component of iron grains with the optical properties 
from Fischera (2004), and the model of Weingartner & Draine 
2001 (see also Li & Draine 2001; Draine & Li 2007). Although 
not unique in terms of reproducing the properties of the galac- 
tic cirrus, these models do nevertheless differ in some predic- 
tions relevant to the modelling of the infrared emission of spiral 



galaxies, particularly with regard to whether conducting parti- 
cles are included as a constituent of the dust model. For exam- 
ple in the model of Fischera & Dopita a substantial part of the 
diffuse 60 //m emission is powered by optical photons absorbed 
by iron grains in equilibrium with the ambient radiation fields 
(see also Chlewicki & Laureijs 1988), whereas in the model 
of Weingartner & Draine (2001) this emission mainly arises 
from carbonaceous grains stochastically heated by UV photons. 
Potentially, therefore, the choice of grain model can influence 
the inferred contributions of young and old stellar populations 
in heating the grains. We are however not aware of any direct 
evidence for the existence of pure metallic grains in the diffuse 
interstellar medium of galaxies. Furthermore, simple considera- 
tions of potential redox reactions indicate that any such grains in 
the diffuse neutral ISM will be susceptible to oxidation and will 
revert to the properties of non-conducting grains (Duley 1980), 
like silicates. Here we substituted the model of Draine & Lee 
(1984) and Laor & Draine (1993) used in Papers I-IV to the 
latest model from Weingartner & Draine (2001) and Draine & 
Li (2007) incorporating a mixture of silicate, graphite and PAH 
molecules. 

We note that the adopted dust model is only used for the cal- 
culation of transfer of radiation through the diffuse interstellar 
medium, where the properties of the grains are reasonably well 
constrained. When including the infrared emission from star- 
forming complexes, which contain dense clouds where grains 
may form ices and other complex compounds, we use a tem- 
plate model SED function which is observationally constrained, 
as described in Sect. 12.81 This sidesteps the uncertainties in the 
optical properties of grains in dense clouds. 

For the parameters of the model we consider Case A and 
Rv = 3.1 from Weingartner & Draine (2001), with the updates 
from Draine & Li (2007) (revised size distribution and optical 
constants for PAHs). The parameters of the model considered 
here are summarised in Table lE.3l 

Following Li & Draine (2001) we allow the relative abun- 
dance of neutral and ionised PAHs to vary with molecule size but 
to be independent of position in the galaxy (or intensity of the ra- 
diation field). The ionisation fraction (<!>,■„„) of PAHs was fixed 
according to the average over the three phases of the diffuse ISM 
given in Fig. 7 of Li & Draine (2001). The PAH-size distribution 
was then multiplied with <t>, on and 1 - <t>, „, respectively, to ob- 
tain the size distribution of ionised and neutral PAH-molecules, 
respectively. 

Since for the dust and PAH emission we accurately take into 
account the stochastic heating of the grains by the ambient ra- 
diation fields (see Sect. 12.6b . we also need to derive the heat ca- 
pacities for our dust grains. The heat capacities for silicate and 
graphite grains needed to derive the temperature distributions are 
summarised in Popescu et al. (2000). The heat capacities for the 
PAH-molecules were calculated using Eq. 15 of Li & Draine 
(2001) and Eq. 5 from Li & Draine (2002). 

Finally, to facilitate the incorporation of future improve- 
ments in the dust model/inclusion of new dust models, we have 
constructed a flexible interface between the code for calculating 
dust emission and the input dust model. 

Because we have modified the dust model we also had to 
redo the calculations of attenuation from Paper III. From the 
perspective of the extinction properties of the two dust models 
(the old version and the new updated version), they are both de- 
signed to reproduce the observed extinction curve of the Milky 
Way. However the discrepancy arises from the change in the rel- 
ative contribution of absorption and scattering to the total extinc- 
tion. This will give rise to a small but non-negligible difference 
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in the overall energy balance (around 10% for NGC 891 - see 
Appendix[B]). In Appendix[B]we briefly describe the attenuation 
calculations obtained with the revised dust model and compare 
with the results from Paper III. 

2.5. Calculation of radiation fields 

In the formulation of our model an explicit calculation of radi- 
ation fields is only done for the diffuse component. There are 
two channels by which the diffuse interstellar medium can be 
illuminated by stellar light, since, as described in Sec. 12.11 our 
model has a diffuse and a clumpy distribution of stars and dust. 
The first channel is through the illumination of the diffusely dis- 
tributed dust by the smoothly distributed population of old stars. 
The second channel is through the escape of radiation from the 
young stars, out of the localised star-forming complexes. 

Thus, the first step in the calculation of radiation fields in the 
diffuse medium is to determine the amplitude and SED of the 
stellar light which escapes from star-forming complexes. 

2.5.1 . The escape fraction of stellar light from the clumpy 
component 

In Sect. [272] we introduced the factor F, which can be physically 
identified with the luminosity-weighted mean fraction of direc- 
tions from the massive stars, averaged over the lifetime of the 
stars, which intersect the birth-cloud. The high opacity assumed 
by our model for the birth clouds bestows a grey body absorption 
characteristic for the locally absorbed radiation in any individual 
star-formation region. There is nevertheless a wavelength depen- 
dence of the probability of absorption of stellar photons because 
stars of different masses survive for different times, such that 
lower mass and less blue stars spend a higher proportion of their 
lifetime radiating when they are further away from their birth- 
clouds, and because of the progressive disruption of the clouds 
by supernovae and stellar winds. To calculate this wavelength 
dependence we follow the analysis of Appendix A of Paper III, 
thus preserving consistency between the model of the attenua- 
tion of light from the young stellar population given in Paper III 
and the model of the infrared emission from star-forming regions 
given here. Denoting the wavelength dependence of the locally 
absorbed radiation by the function Fx, the spectrum and ampli- 
tude of the starlight from the young stellar population in the thin 
disk with spectral luminosity density Lf sk which is intercepted 
by a parent cloud is thus given by: 

jjocal _ x jjdisk (21) 

where 

F,\ — F x fx (22) 

and fx is tabulated in Table A. 1 of Paper III. As described in 
Paper III, the tabulated values of fx correspond to a secular de- 
crease of the solid angle subtended by the parent clouds on a 
timescale of 3 x 10 7 yr. For completeness we remind the reader 
that ]Jf sk is defined in terms of the free parameter of the model, 
SFR, in Eq.[[9](see Sect. l2331 l.FI 

The amplitude and spectrum of the light from the young stel- 
lar population which illuminates the diffuse dust, is given by 

L tdisk,diff _ l _ jlacal (23) 

8 We note that L'f sk and Vf" k denote the same physical quantities, ex- 
pressed respectively as spectral densities in wavelength and frequency. 



Since Fx must satisfy Fx < 1.0 at all wavelengths, our model 
has an intrinsic constraint on the dumpiness factor of 

F < 1.0/max(fx) = 0.61 (24) 

Physically, this constraint corresponds to the requirement 
that the "porosity factor" po of the birthcloud, used in Eq. A2, 
A12 and A13 of Paper III to denote the fragmentation of the 
cloud, is always less than unity. This limit denotes a complete 
lack of fragmentation of the cloud due to the action of super- 
novae and stellar winds, such that the probability of escape of 
photons into the diffuse ISM is determined only by the migra- 
tion of stars away from their birthclouds. 

2.5.2. Combining radiation fields from the young and old 
stellar populations 

Since all the diffuse stellar components are attenuated by the 
same distribution of dust, the diffuse radiation fields seen by each 
grain can be considered to be a sum of the radiation fields pro- 
duced by each stellar component. Thus we calculated the diffuse 
radiation fields separately for the three main stellar components 
of our model, the young stellar disk, the old stellar disk, and the 
bulge. Then the individual diffuse radiation fields can be com- 
bined to produce the total diffuse radiation fields of our desired 
model galaxy. This is the same concept as that used in Paper III 
for the calculation of the attenuation of stellar light. 

To achieve this we first ran the radiative transfer calcula- 
tions for each stellar component, each wavelength and each 
given optical depth r B , using the intrinsic spectral luminosity 
densities L okl ., and L youn ?,, as tabulated in Table IE. 21 This re- 

v, unit v, unit ' 1 1 

suited in the calculation of the unit energy densities for the disk 
uff nit (R, z, T^j), for the thin disk uf sk nit {R, z, r£) and for the bulge 
bulge ,p f n j. ^ unit stellar luminosities (i.e. L old ., for the 

A, unit - ' ' S y v y, unit 

disk or the bulge and l/ am ? for the thin disk). 

& v, unit ' 

Then the unit radiation fields were scaled according to the lu- 
minosity of each diffuse stellar component, according to the pa- 
rameters of the model. For the old stellar disk (only diffuse) the 
resulting radiation fields uf k (R, z, t b , old) were obtained from: 

uf k {R, z, t b , old) = old x uf k nit (R, z, r f B ) (25) 

where old is the normalised luminosity of the old stellar 
population, defined by Eq. [T2] For the young stellar disk 
uf sk (R, z, t£, S FR, F) is derived from: 

uf sk (R, z, r£, S FR, F) = f FR x(l-Ff A )x 
lM yr _1 

xiO**.^) (26) 

where S FR is related to the normalised luminosity of the young 
stellar populations in Eq. [16] For the bulge the radiation fields 
are given by: 

u A ulge (R, z, r f B , old, B/D) = old x B/D x uf** t (R, z, r£) 

(27) 

The total diffuse radiation fields 

ux(R,z, t^ b ,S FR, F, old, B/D) are given by: 

u A (R, z, t^ b ,SFR, F, old, B/D) = uf k (R, z, t{, old) + 
+uf sk (R, z, t^ b ,SFR, F) + u"' ge (R, z, t{, old, B/D) (28) 
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This formulation provides an exact treatment of the spatial 
variation of the colour of the radiation field in galaxies, which is 
a crucial factor in shaping the infrared SEDs. 

In principle this formalism implies that we would need to 
create a library of radiation fields calculated for each possible 
value of the F factor, implying a different wavelength depen- 
dence (1 - Ffx) of the fraction of photons escaping from the 
clumpy component into the diffuse one, making the problem 
quite intractable. In practice, however, it is possible to simplify 
the problem if we make the assumption that the wavelength de- 
pendence (1 - Ffx) is fixed to the (1 - F ca ifx) corresponding to 
a fiducial calibration value F ca \ = 0.35, which we consider to 
be a typical value for spiral galaxies in the local universe. This 
intermediate value of F ca \ (theoretical values of F can vary from 
to 0.61; Paper III) is consistent with the observation that 60- 
70% of the dust emission from spiral galaxies emanates from 
the diffuse component (e.g. Sauvage et al. 2006). Furthermore, 
detailed analysis of face-on spiral galaxies with resolved star- 
formation regions in the Spitzer 24 /mi and the GALEX FUV 
bands (e.g. Calzetti et al. 2005) have shown an almost one- 
to-one correspondence between the direct and re-radiated light 
coming from these complexes, albeit with a large variance in 
the ratio of the two emissions. Thus there were essentially no 
Spitzer sources with no FUV counterparts (meaning that F is 
unlikely to take the value 0.61), but also no FUV sources with 
no Spitzer counterparts (meaning that F is unlikely to take the 
value 0). We note that the value of F ca i is only used to determine 
the exact wavelength dependence of that part of the light from 
the young stellar population which escapes the star-formation 
complexes and actually illuminates the diffuse dust. Although 
the computed colour of this light will be inaccurate if F differs 
from F ca i, corresponding inaccuracies in the predictions for the 
dust emission SEDs are unlikely to be more than a few percent, 
as demonstrated in Sect. 14.31 where we quantify the effect on 
the dust emission of changing the spectral shape of the stellar 
emissivity. We note that, for the same reason, when dereddening 
the observed UV/optical SED, the wavelength dependence of the 
component of attenuation from the clumpy component can use 
the exact solution for F found in the optimisation. 

The effective reddening law of the continuum photons escap- 
ing the star-formation regions corresponding to F ca i is given in 
Table EH 

The use of the calibration factor F ca i means that in practice 
all the equations that describe the illumination of the diffuse dust 
by the young stellar disk need to be rescaled to accommodate 
different values of F than those used in the calibration. For this 
we define a correction factor for the diffuse component corn(F): 



in which case Eq.l2T1becomes: 



corr d (F) = (1 - F) X 



j young 
unit, UV 



(29) 



In this case Eq.l26*lbecomes: 



uf s \R, z, r f B ,SFR, F) = f FR (1 - F ca if a) cor/(F) x 
D 1 Mo vr _l 



<*O*,z.*£)(30) 



In an analog way we define a correction factor for the localised 
component corr'(F): 



corr'(F) -Fx 



jjoung 
unit, UV 



JUV A, unit tin j a 



(31) 



r local 



SFR 
1 M yr 



-F cal f A corr'(F)xL 



young 
A, unit 



(32) 



2.5.3. The radiative transfer calculations 



One of the key elements of any self-consistent model for the 
SEDs of galaxies is the incorporation of a radiative transfer (RT) 
code that can be used to determine the radiation fields in galax- 
ies and their appearance in the UV/optical bands. Four meth- 
ods have been proposed so far in RT models, though only two 
of them have been broadly used in the context of galaxy mod- 
elling, namely ray-tracing methods (Kylafis & Bahcall 1987, 
Silva et al. 1998, Semionov & Vansevicius 2005a, Semionov 
& Vansevicius 2006) and Monte-Carlo (MC) techniques (Witt 
et al. 1992, Bianchi et al. 1996, de Jong 1996, Gordon et al. 
2001, Baes & Dejonghe 2002, Baes et al. 2003, Baes et al. 2005, 
Jonsson 2006, Baes 2008, Bianchi 2008, Chakrabarti & Whitney 
2OO8).0 

Since our model makes explicit calculations of radiation 
fields in galaxies we use an ray-tracing code, which is better 
suited to this type of problem. As in Paper I we use the Kylafis & 
Bahcall (1987) code, updated and revised according to the spe- 
cific needs of our modelling technique. 

The code of Kylafis & Bahcall (1987) was originally de- 
veloped to model the surface brightness distribution of edge-on 
galaxies with cylindrical symmetry. This utilises a ray-tracing 
method that exactly calculates the direct and first order scattered- 
light and uses the method of scattered intensities to approximate 
the higher order scattered light (see review of Kylafis & Xilouris 
2005). Several modifications were incorporated for the present 
application, the most significant of which we briefly describe 
here. 

Firstly, the code used was generalised and optimised for ap- 
plication to any cylindrically symmetric distribution of emissiv- 
ity and opacity, such as the thin disk representing the young stel- 
lar population and associated dust in our standard calculations, 
and the annular representation of this disk that we use in the 
tests made in Sect. 14.21 to evaluate the second disk approxima- 
tion to spiral structure. To achieve this, we incorporated an in- 
ray sampling scheme which automatically adapts to variations 
in emissivity and opacity encountered along the rays. Further, in 
regions of low emissivity, where the method of scattered intensi- 
ties becomes less accurate, the maximum order of the scattered 
light had to be reduced. The reduction was done progressively 
accordingly to the emissivity encountered, such that only the 
fully accurate first order scattering term was utilised in the re- 
gions with the smallest emissivity. This is a good approximation 
for the geometries of dust and stars considered here, in which 
low emissivities are never encountered in regions of high opaci- 
ties. Lastly, for application to the calculation of radiation energy 
densities the code was further modified to calculate the vector 
radiation field sampled in each direction of a predefined set of 
216 directions. 

To provide a more frequent sampling of the inner part of the 
galaxy the spatial positions were arranged on an irregular grid 
with sampling intervals in radial direction ranging from 50 pc in 



9 There are also RTs initially developed for use in other astrophysical 
contexts (Mattila 1970, Witt 1977, Yusef-Zadeh 1984, Wood et al. 1996, 
Lucy 1999, Wolf et al. 1999, Bjorkman & Wood 2001, Steinacker et al. 
2003, Wolf 2003, Juvela 2005, Steinacker et al. 2006, Fritz et al. 2006). 
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the centre to up to 2 kpc in the outer disk and with sampling inter- 
vals in vertical direction ranging from 50 pc in the plane to up to 
500 pc in the outer halo. This scheme was necessary to properly 
sample the radiation in the more optically thick central regions, 
where, as can be seen from Fig. [4] the dust associated with the 
thin stellar disk in conjunction with the strong positional depen- 
dence of the emissivities of the young stellar population (in z) 
and the inner emissivity cusp of the bulge (in both R and z) can 
lead to quite small scale structures even in the diffuse radiation 
fields. We made some further tests with higher sampling, rang- 
ing from 25 pc to up to 1 kpc in radial direction and from 25 pc to 
up 250 pc in vertical direction. This made a less than 1% differ- 
ence in the results with a large increase in the computation time. 
We have therefore concluded that the optimum sampling of the 
radiation fields is the AR - 50 pc - 2 kpc and Az = 50 pc - 500 pc 
scheme. 

The database of integrated SEDs presented in this paper 
does not incorporate the effect of the reabsorption of infrared 
light emitted by grains in the galaxy, as this process has only a 
marginal effect on the integrated dust emission, and would en- 
tail a significant increase in complexity and decrease in speed 
of the calculations. This is a good approximation for the present 
problem. For example, for the most opaque disks we consider 
(r^ = 8), only the ~ 3 /j.m dust emission shows a significant 
depression in the integrated emission. However we do calculate 
reabsorption of infrared light when we produce maps of galaxies 
seen at higher inclination, and this process has been incorporated 
in the calculation of the profiles presented in Sect. 14.11 

The radiation fields illuminating the dust in the diffuse disk 
are calculated taking into account anisotropic scattering, using 
the albedo and anisotropy of the scattering phase function as 
defined in Sect. [A] and the geometry of the dust as defined by 
Eqs. [4] and [6] in Sect. [2J] At each UV/optical wavelength^ of 
Table IBT21 and each optical depth r B a separate calculation was 
made for each of the three stellar components, with smoothly 
distributed emissivities given by Eq.[T]for the old stellar popula- 
tion in the disk and in the bulge, and Eq.[5]for the young stellar 
population in the thin disk. This adequately samples the wave- 
length variation in the optical constants of grains; tests done on 
this indicate uncertainties at the level of a percent. Concerning 
the corresponding error due to the finite sampling of the stellar 
SEDs, the reader is referred to Sect. 14.31 in which we show that 
the exact colours of either young or old component of the stellar 
SEDs only affect the resulting dust emission SEDs at the level 
of a percent. Overall, therefore, we consider that uncertainties 
introduced by our finite wavelength sampling on the predicted 
dust emission SEDs are minor. 

When calculating the heating of the grains, and therefore of 
the wavelength-integrated energy absorbed, the radiation fields 
were interpolated on a fine wavelength grid containing 500 data 
points logarithmically sampled between 912 and 50000 A. Tests 
done with 1000 data points showed negligible differences in the 
results, therefore we decided to fix the wavelength sampling to 
500. Energy balance tests done within a common wavelength 



sampling indicate that our models conserve energy to within bet- 
ter than 1 percent. 

2.6. Calculation of grain temperature distributions 

As described in detail in Popescu et al. (2000) our calculation 
incorporates an explicit treatment of the temperature fluctua- 
tions undergone by small grains whose cooling timescales are 
shorter than the typical time interval between impacts of UV 
photons. This so-called "stochastic heating" process determines 
the amplitude and colour of the bulk of the diffuse emission from 
most spiral galaxies in the shortest infrared wavelength bands 
(< 100 /im). Thus, since for each position in the galaxy the en- 
ergy density of the radiation fields heating the grains has been 
derived (see Sect. 12.5b , we can calculate the probability distribu- 
tion of temperature P,(a, T) for each grain size a and composi- 
tion i by equating the rate of absorption of energy to the rate of 
emission of energy : 



7i a 2 c J Q abSy , 
: J B A (a, 



j(a,A)u^dA = 4 
x | H ; (a,T)P i (a,T)dTdA 



n(na 2 ) J" 



Qabs,i(a,A) x 



(33) 



where c is the speed of light and Bx is the Planck function. 

The probability distribution of temperature is calculated fol- 
lowing Voit (1991), which combines the numerical integration 
method of Guhathakurta & Draine (1989) and a stepwise analyt- 
ical approximation. We sample the temperature distribution in 
61 logarithmically spaced points. 

2. 7. Calculation of infrared emissivities and spatially 
integrated SEDs in the diffuse component 

Once the temperature distribution of a grain is known, the bright- 
ness 1^ i of a grain of radius a and composition i is then given by: 



h,i(a) = Qabsj(a,A) B A (a,T)Pi(a,T)dT 



(34) 



where I^.i is in units of [Wm^sr'A -1 ]. 

Once the infrared brightnesses , are calculated for each 
grain size and composition, we can then integrate over the grain 
size distribution n{a) to obtain the infrared emissivity per H atom 
of grains of a given composition i, .: 



na n(a) h,i(a) da 



(35) 



where f? . is in units of [W sr 'A 'H 1 ]. 

Finally, the total infrared emissivity per H atom, is ob- 
tained by summing over the grain composition i 



10 The radiation transfer calculations were performed at each of the 
wavelengths given in Table lE,2l with the exception of the longest wave- 
length at 50000 A, where the radiation fields could be extrapolated suf- 
ficiently accurately from the radiation fields at 22000 A assuming a 
Rayleigh- Jeans law. The sampling points for the calculations are the 
same as those used for the calculation of attenuation in Paper III, except 
for the addition of points at 1500 & 3650 A to improve the definition of 
the stellar emissivity SEDs. 



./!/ 



(36) 



In Fig. [2] we show the calculation for the infrared emis- 
sivity of grains heated by the local interstellar radiation fields 
(LIPsF), with the energy densities of the radiation fields taken 
from Mathis, Metzger & Panagia (1983). The resulting infrared 
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Fig. 2. The dust and PAH emission SED for the diffuse ISM at 
high Galactic latitude (solid line) calculated using our model. 
Also plotted here are the COBE data (symbols with error bars). 
These are given as an average of the data from the North eclip- 
tic pole field and the Lockman Hole field (Arendt et al. 1998) 
and are further colour corrected. The contributions of the differ- 
ent dust compositions to the total model SED are as follows: Si 
(dotted line), Gra (dashed line), PAH (dashed-dotted line). 

SED is similar to that presented in Fig. 12 of Draine & Li (2007) 
and indicates that our method for calculating stochastic heating 
of grains gives similar results to the method of Draine & Li. 
This can be taken as a benchmark test for the calculation of in- 
frared emission. In Fig. [2] we also plotted the colour-corrected 
DIRBE/COBE data for the high Galactic latitude. As expected, 
there is a good match between the model and the observations. 
In Fig.|2]we also show the predicted contribution of the different 
dust compositions to the total infrared emissivity. 

Once the total infrared emissivities per H atom have been 
calculated, these can be then scaled to the volume density of 
dust grains at each position in the galaxy to obtain the volume 
luminosity density of the first and second dust disks, Lff (R, z) 

and L'f'j'' (R,z), as a function of position: 



Lff ust (R,z) = 



An 



t-> 't ■ 
1 g 1 ratio 



C ext (B) 2zf k (l+T ra , !0 ) 



x #(*,z)x 



x exp 



R 

kdisk 
"d 



\Z\ 



(37) 



j tdisk 
^A,dust 



(R,z) = 



4tt 



CeAB) 2 Z f sk (l+T r , 



x exp 



x /?(/?, z)x 



R 

U tdisk 
"d 



Izl 



(38) 



where L dtsl 



and V 



tdisk 



are in units of [Wm 3 A '] and 



X,dust A,dust 

Cext(B) = 6.38 x 10 _22 cm 2 H _] is the extinction cross-section 
in the B band (at 4430 A) for the dust model considered here. 
The spatially integrated SEDs (the spectral luminosity density) 
of the diffuse component of the model galaxy (with the fixed size 



diff, model 



/i disk ,), L\ , 

s, ref J ^ A, dust 



is then: 



j diff, model 
"'A, dust 



, j tdisk 
'A,dust AAust 



)RdRdz 



(39) 



We checked the fundamental energy balance in the calcu- 
lations, by comparing the energy absorbed from the radiation 
fields with the energy emitted in the infrared, obtaining consis- 
tency to within better than 1%. 

2.8. The infrared emission from the star-forming regions 

In our model we take the total energy absorbed by the birth- 
clouds to be 



r local 
-'abs 



I rhc, 
J A 



r tdisk 



wn-uv ^ion-uv 



(40) 



V*£l m is defined in Eq.QUand 



where jj? cal is defined in Eq. 

fion-uv — 0.3. 

That is, as justified in Sect. 12.21 we take the heating of the 
grains in the birthclouds to be dominated by photons from the 
stellar progeny, and neglect any external contribution from the 
ambient radiation fields in the galaxy. 

To determine the SED of the reradiated light from dust grains 
and PAH molecules in the birthclouds we use a fixed spectral 
template of a photodissociation region (PDR): 



j local _ j local j PDR 
^A,dttsl ~ ^abs ^A 



(41) 



where L P X DR is the template function for a photodissociation re- 
gion, normalised such that 



L p x DR dX 



1.0 



(42) 



We use the model of Groves et al. (2008) (see also Dopita 
et al. 2005), fitted to broad-band measurements of dust and PAH 
emission from galactic star-formation regions, to provide an em- 
pirically constrained prediction for the detailed spectral form of 
the template L PDR . Although this model was primarily developed 
for use in predicting the SED of starburst galaxies, its fundamen- 
tal constituent is a prediction of the PAH/dust and nebular line 
emission from the HII region and PDR components of individ- 
ual SF regions, and so, with suitable choice for the values of 
the model parameters, is also directly applicable to the predic- 
tion of the PDR emission from SF regions in spiral galaxies. For 
the spectral template used here we only consider the PDR com- 
ponent calculated by Groves et al., ignoring the emission com- 
ponent from the HII region (as justified below). We also only 
include the PAH and dust emission components in the spectral 
template, not including the free-free and line emission from the 
gas phase. 

The prime physical motivation for the use of the model of 
Groves et al. is that it self-consistently calculates the effect on 
the emergent dust and PAH emission SED of the dynamical evo- 
lution of the emitting regions. Specifically, the model considers 
a spherically symmetric, fully enclosed mass-loss bubble, driven 
by the mechanical energy input through winds and supernovae of 
a star cluster as a function of the external density (parameterised 
in terms of the pressure of the ambient ISM) and a "compactness 
parameter", C, which scales according to the mean luminosity- 
weighted flux of photons onto the PDR. The emission from PAH 
and dust in the PDR is then self-consistently determined through 
a separate radiation transfer calculation (Groves et al. 2008), as 
a function of the grain column density in the PDR (which in 
the formulation of Groves et al. is actually prescribed through a 
combination of the gas column in the PDR and the metallicity). 
Being based on a dynamical model, the calculations of Groves 
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Table 1. The references used for collecting the observed flux 
densities of the 57 SF regions plotted in Fig. [3] 



Bains et al. (2006) 

Buckley & Ward-Thompson (1996) 

Conti & Crowther (2004) 

Gordon (1987) 

Gordon & Jewell (1987) 

Hill et al. (2005) 

Mooney et al. (1995) 

Moore et al. (2007) 

Sievers et al. (1991) 

Thompson et al. (2006) 

Ward-Thompson & Robson (1990) 



et al. explicitly take into account the expected variation in colour 
and amplitude of the incident photon flux on the PDR due to the 
increase in radius and evolution of the stellar population, which 
acts to broaden the SED of the dust emission from the large en- 
semble of SF regions at different evolutionary stage expected in 
a spiral galaxy. The model also considers the photodestruction 
of PAH, strongly influencing the relative amplitudes of the PAH 
and far-IR/submm emission components from the PDR. 

To empirically verify these predictions of the model, and de- 
termine the best model parameter to determine the PAH and dust 
emission template from an ensemble of localised SF regions in 
spiral galaxies, we compared the model predictions with ob- 
served data (see Table [TJ for the radio-selected sample of star- 
formation regions in the Milky Way of Conti & Crowther (2004). 
This sample is distributed throughout the disk of the Milky Way, 
so is likely to be representative of the population of SF regions 
in the galaxy. In Fig. [3] we plot the mid-IR/far-IR/submm SEDs 
of 57 galactic star-forming regions from Conti & Crowther, nor- 
malised to the 100 yum flux density. Overlaid on this is the best fit 
model for the PDR emission from Groves et al., found by search- 
ing the parameter space in the compactness factor C and the dust 
column density in the PDR (adjusted by varying the product of 
metallicity and gas column density). The other major parameter, 
pressure, was kept fixed as the PAH and dust emission is com- 
pletely insensitive to this parameter which only affects the emis- 
sion lines (see Fig. 4 of Groves et al.), which are not required in 
our template. In comparing the model with the data we did how- 
ever take into account the free-free emission component from the 
PDRs, since this emission component will also be represented in 
the broad-band continuum data of the Milky Way star formation 
regions. The free-free emission slightly raises the level of the 
emission near 1 mm above the prediction of the dust emission 
template, as shown by the divergence of the dotted line, which 
includes the free-free component from the full line in Fig. [3] 

The best fitting PDR model is for compactness parameter 
log(C) = 6.5, metallicity 1.0, and hydrogen column density 
log(N) = 22.0. Although there is a large spread in observed 
colours, indicating the wide range of evolutionary states and en- 
vironment of individual SF regions, overall this model repro- 
duces the mean colours of the population of sources over the 
whole spectral range covered by the data. This spectral range 
extends from 12yum (including the PAH emission), through 25, 
60 and lOO/vm (where emission is dominated by large grains in 
thermal equilibrium with the intense local radiation fields), and 
into the submm/near-mm range (where the emission is domi- 
nated by dust which is cold due to the self-shielding of grains 
in the PDR. This latter cold dust emission component is difficult 
to predict theoretically, depending on the otherwise poorly con- 
strained dust column density, so is fundamentally an empirically 
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Fig. 3. MSX and IRAS mid-IR and far-IR integrated photometry 
for 57 SF regions from Conti & Crowther (2004), supplemented 
by all published sub-mm and/or near-mm data for which inte- 
grated flux measurements covering the full angular extent of the 
source are available (18 sources); (see Table [TJ. Measurements 
are normalised to the 100//m flux density. The solid line is the 
adopted model for the PDR dust, PAH emission from Groves et 
al. (2008), having compactness parameter log(C) = 6.5, metal- 
licity 1.0, and PDR column density log(N) = 22.0. The dashed- 
line, which deviates slightly from the solid line in the submm 
range, is the same but with the free-free emission included. 

determined quantity. Likewise, the relatively high value of the 
compactness factor C is needed to fit the rather warm observed 
far-IR/mid-IR colours. As directly confirmed through high res- 
olution imaging of the sources this indicates a close proxim- 
ity (typically parsec scales) of the exciting star clusters to the 
PDRs, most of which are associated with rather dense fragments 
of molecular material which will be only slowly pushed out from 
the cluster. It is interesting to note that the fit of the PDR model 
to the data leaves no room for a warmer dust emission compo- 
nent from an enclosed HII region, suggesting that the HII regions 
typically extend to larger spatial scales. On this picture, the cloud 
fragments carrying the PDR emission are embedded in a diffuse 
ionised medium which blends in with the diffuse emission from 
the second dust disk. 

Although, as we argue above, the physical values for the 
compactness parameter C, metallicity and hydrogen column 
density used for the fit seem plausible, we note that there may 
be room for some systematic uncertainties in these parameters 
if the optical properties of the grains differ from the model of 
Weingartner & Draine used by Groves et al.. The most important 
point however in the context of our use of this template SED is 
that it gives a good fit to a representative sample of star formation 
regions in our Milky Way. So even if the physical interpretation 
of the shape of the SED may be still somewhat imprecise, the 
template itself will be the correct one. In conclusion, we believe 
that the PDR model of Groves et al. that best fits the galactic star- 
forming regions can be taken as a template SED for the clumpy 
component of our model. 

2.9. The free parameters for the calculation of the infrared 
SEDs 

The spatially integrated SED for the dust emission of a galaxy is 
given by: 

L™% l st (T f B ,SFR,F,old,BID) = 
L fju7 de1 ^ s FR > F > ° ld > B / D ) + *%LV FR > F ) ( 43 ) 



Cristina C. Popescu et al.: Modelling the spectral energy distribution of galaxies. 



15 




Fig. 4. Examples of calculated radial profiles of radiation fields 
in the plane of the disk (z — pc), for the best fit model of our 
prototype galaxy NGC 891. The different profiles are for differ- 
ent wavelengths. 

where L d //;™ del is defined by Eq.|39]and L'^JSFR, F) is de- 
fined by Eq.|4TJ From Eq.|43]one can see that the free parameters 
of our model are, to recap: 

- the central face-on opacity in the B-band (Eq. [TT1 
Sect. ED, 

- the dumpiness factor F (Sect. I2.2K 

- the star-formation rate SFR (Sect. l2~33T >. 

- the normalised luminosity of the old stellar population old 
(Sect. 12311 and 

- the bulge-to-disk ratio B/D (Sect. I2T21 . 

In cases where detailed modelling of the optical/NIR images is 
also available, e.g. the study of NGC 891, the old and B/D pa- 
rameters are independently constrained from the optimisation of 
the optical images, and thus only 3 free parameters are needed 
for the model: t b , SFR, F. 

We note that the free parameters r^, F, and B/D are also 
free parameters for the attenuation of the UV/optical stellar light 
given in Paper III (the fourth and final free parameter that affects 
the attenuation is the inclination i of the disk). Therefore the 
dust emission SEDs predicted here can be used in conjunction 
with the predictions for attenuation given in this paper for a self- 
consistent modelling of the UV/optical-IR/submm SEDs. 

The parameter t b can be related to the total dust mass using 
the equation: 

M dusl = yh](B)T f B (44) 

where y is a constant related to the geometry of the distribution 
of dust in galaxies and to the dust model. For our model with 
h]{B) in parsec, y = 0.9912M G pc~ 2 . 

3. Illustration of the model on NGC 891 

We illustrate our model, including description of the interme- 
diate steps in the calculation, using NGC 891. Because of the 
changes to our model, it is important to see if we can now fit the 
whole SED of NGC 891, including the MIR emission dominated 
by the PAH spectral features. 



Fig. 5. Examples of calculated vertical profiles of radiation fields 
in the centre of the disk (R = Opc), for the best fit model of 
our prototype galaxy NGC 891. The different profiles are for 
different wavelengths. 

As described in Sect. 12.91 since the luminosity of the old 
stellar populations has been derived from the optimisation of the 
optical/NIR images, we only need 3 free parameters to fit the 
dust emission SED of NGC 891: r^, SFR, F. For this we ran 
the calculations for trial combinations of these parameters. We 
note that the normalised luminosity of the old stellar population 
derived for NGC 891 is old = 0.792. Throughout this section all 
the examples shown are for the best fit modeQ for NGC 891, 
which is for t£ = 3.5, SFR = 2.88 M /yr and F = 0.41. 

3.1. The radiation fields 

The first step in the calculation (see Sect. I2.5l l is the derivation 
of the radiation fields illuminating the diffuse dust, which we 
illustrate in Figs. [4] and [5] through examples of radial and verti- 
cal profiles of energy densities of total radiation fields (not 
to be confused with profiles of stellar emissivities). The sharp 
rise in the inner parts of the radial profiles in the K and B band 
is produced by the dominance of the radiation coming from the 
bulge. We note here the large variation in the colour of the ra- 
diation fields with position, in particular in the radial direction, 
which, as shown in the next figures, will introduce large differ- 
ences in the shape of the FIR SEDs. Thus, models that assume 
radiation fields with the fixed colour of the local interstellar radi- 
ation fields (LIRF) are likely to introduce systematic uncertain- 
ties in the predictions for the dust emission SEDs. 

3.2. The temperature distributions 

The next step in the calculations (see Sect. I2.6l l is the derivation 
of the probability distributions of dust temperature. In Fig. [6] we 
show examples of temperature distributions for grains of differ- 
ent sizes and compositions placed at different positions within 
the galaxy. The sizes plotted for each composition reflect the 
range of sizes given by the dust model considered in this paper. 
For example sizes smaller than 10A are only considered for PAH 

11 We quote the best fit S FR to two decimal places because we use 
S FR as a proxy for the luminosity of the young stellar population (see 
Eq.[[6j. 
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molecules. Also, according to the dust model, the biggest PAHs 
are those of 100 A, while graphites have only sizes larger than 
100A. As expected, small grains exhibit broad probability distri- 
butions. For example PAH molecules are only heated stochasti- 
cally, as they have small sizes. Larger grains have narrower prob- 
ability distributions, eventually tending towards a delta function, 
when they start to emit at equilibrium temperature. 

The figure also illustrates the dependence of the temperature 
distributions with position in the galaxy, due to the spatial varia- 
tion of the intensity and colour of the radiation fields. For exam- 
ple a graphite grain of size 0.0316 yum will exhibit a delta func- 
tion temperature distribution if placed in the centre of the galaxy 
(first and third row from the top, middle panels, blue curves in 
Fig. [5}, where the radiation fields have higher energy densities 
and redder colours, than if placed in the outskirts of the galaxy 
(second and 4th rows from the top, middle panels). Thus, in the 
centre of the galaxy the 0.0316/im graphite grain will emit at 
equilibrium temperature, while the same grain will start to emit 
stochastically at/? = 15000 pc. 

Finally Fig. [6] shows the dependence of temperature on com- 
position. For example a 0.0316 /mi grain placed in the centre of 
NGC 891, will exhibit a lower temperature if it had a silicate 
composition (left upper panel) than if it had a graphite compo- 
sition (middle upper panel), due to the differences in the optical 
properties of the grains. 

3.3. The infrared emissivities 

In Fig.[7jwe show examples of calculated infrared brightnesses 
(from Eq. [34l l for grains of different sizes and compositions 
placed at different positions within the diffuse ISM of our proto- 
type galaxy NGC 891. For the calculations we used the proba- 
bility distributions of temperature shown in Fig. [6] We also show 
the same positions and sizes as those from Fig. [6] Following the 
trends in the temperature distributions, one can see the depen- 
dence of the infrared brightnesses with position in the galaxy, 
due to the variation in the intensity and colour of the radiation 
fields. One should note that here we show the absolute bright- 
ness of each grain, meaning that the emission is not weighted to 
take into account the abundance of grains according to size or 
composition; it simply indicates the response of each dust grain 
to the radiation fields. 

By comparing the pairs of positions, first and third row from 
the top on one hand, and second and 4th row on the other hand, 
we see that overall the SEDs show a stronger radial variation 
than a vertical variation and that this trend is independent of 
grain composition and grain size. This is primarily a result of the 
fact that the radiation fields show a stronger variation in colour 
in the radial direction than in the vertical direction (see again 
Fig. |5] and |4]i, as described in Sect. 13.11 which, in turn, is a direct 
consequence of the finite disk plus bulge description of our 
geometrical distributions of stars and dust. We will discuss here 
both the change in the peak of the SEDs as well as the change 
in the overall amplitude of the infrared brightness. As expected, 
for a given grain size and composition the peak of the SEDs 
shifts towards longer wavelengths and its amplitude decreases 
in weaker radiation fields, especially at large galactic radii. 
However the wavelength shift of this peak and its amplitude as 
a function of grain size have a more complex behaviour. 

i) The wavelength variation of the peak of the infrared bright- 
nesses with grain size for a fixed galaxian position ( radiation 
field). 



For the stochastically heated grains there is a strong shift in 
the peak of the SED towards longer wavelengths with increas- 
ing grain size. This is seen for the small grains (Si and PAH; 
the green curves) placed at large galaxian radii (second and 4th 
row), which all exhibit strong stochastic heating (as also seen 
from the broad probability distribution of temperature in Fig. [6] 
corresponding panels). Since most of the energy is radiated at 
the highest temperature side of the probability distributions, and 
since the increase in the grain size will decrease the width of the 
probability distributions, this means that grains with bigger sizes 
will reach systematically lower maximum temperatures in the 
probability distributions, thus radiating at longer wavelengths. 

By contrast, grains in equilibrium temperature with the ra- 
diation fields will show a small shift in the SED peak with in- 
creasing grain size, and in the opposite direction, namely to- 
wards shorter wavelengths (see the blue curves from the first 
and third row from the top, left column of Fig. [7]). According to 
Eq. [34] the peak of the infrared brightness will be determined 
by the wavelength dependence of Q a i, s in the far-infrared and by 
the equilibrium temperature at which they radiate. For a given 
radiation field u ra d, the equilibrium temperature depends only 
on the Q a b s (see Eq.[33l). Since bigger grains are more opaque at 
shorter wavelengths, tending to a black-body case, they absorb 
more efficiently, and therefore their equilibrium temperature is 
higher, providing the Q„i, s has the same wavelength dependence 
in the far-infrared, independent of grain size. This is indeed the 
case for silicate grains. The graphite grains however have some 
variation in the wavelength dependence in the far-infrared with 
grain size, which produce the non-monotonic shift in the peak 
of their infrared brightness (see first and third row from the top, 
middle column from Fig. [7J where the blue curves intersect). 

A special case is represented by the smallest PAH molecules 
(the red curves on the right column of Fig. [7). Their infrared 
brightnesses are dominated by the emission bands due to 
vibrational transitions, and these features occur at the same 
wavelengths, independent of the molecule size. This is because 
the vibrations seen in the mid-infrared correspond to the 
fundamental stretching or bending modes of the C-H and C-C 
bonds, and do not involve the molecule as a whole. 

ii) The wavelength variation of the peak of the infrared bright- 
nesses with grain size for a variable galaxian position ( radiation 
fields). 

The shift in the peak of the infrared brightness towards longer 
wavelengths with increasing grain size for stochastically heated 
grains (described at z) is a strong function of galaxian position. 
By comparing the first and 3rd rows from the top with the sec- 
ond and the 4th rows for the small Si and PAH grains (green 
curves, Fig. [7J we see that the shift becomes less pronounced 
for grains at small galaxian radii, where the radiation fields are 
stronger and redder. For the largest PAH molecules (or corre- 
sponding silicate grains) there is almost no shift, their peak re- 
maining constant with wavelength. This shows that these grains, 
despite being small, are in a transition phase towards equilibrium 
temperature, as also proven by their narrower temperature distri- 
bution from the corresponding panels in Fig. |7j So even PAH 
molecules can start to emit closer to equilibrium temperature if 
placed in the centre of the galaxy. 

Conversely, the described shift of the peak of the SEDs to- 
wards shorter wavelength with increasing grain size for grains 
heated at equilibrium temperature is also a strong function of 
galaxian position. By comparing the first and 3rd rows from the 
top with the second and the 4th rows for the big silicate and 
graphite grains (blue curves) we see that the shift can be reversed 
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Fig. 6. Temperature distributions for dust grains of different sizes (plotted as different curves in each panel) and various composition: 
Si (left panels), Gra (middle panels) and PAH + (right panels), heated by the diffuse radiation fields calculated for the best fit model 
of our prototype galaxy NGC 891. Temperature distributions for PAH are not plotted in this figure. The colour codding is as 
follows: red is for grains with radius a < 0.001 fxm (0.00035, 0.00040, 0.00050, 0.00063 and 0.00100 jum), green is for grains with 
0.001 < a < 0.01 (0.00158, 0.00251, 0.00398, 0.00631, 0.OlOOO^m) and blue is for grains with a > 0.01 (0.0316,0.10000,0.31623, 
0.7943 /im). The biggest grains have delta function distributions, since they emit at equilibrium temperature. Going from the top to 
the bottom panels the calculations are done for different positions in the model galaxy: R = Opc, z = Opc; R = 15000 pc, z = Opc; 
R = Opc, z = 600pc; and/? = 15000pc, z = 600pc. 



in the weaker and bluer radiation fields at large galactic radii. 
This means that even big grains can start to emit stochastically 
in the outer regions of galaxies. Indeed, by looking at the cor- 
responding panels from Fig. [6] we see that the big grains have a 
non-negligible width in the probability distribution of tempera- 
ture for the outer disks. 

To conclude, the shift of the peak of the infrared brightness 
as a function of grain size strongly depends on the temperature 
of the grains, and therefore on the stochastic or non-stochastic 
nature of the heating mechanism. Since the heating mechanism 
depends both on grain size and on the intensity and colour of the 
radiation fields, it is clear that the shift cannot be described in 



terms of grain size only. This also shows that models that have 
a fixed grain size for the transition between the main heating 
mechanisms of dust, irrespective of the radiation fields, will lead 
to systematic spurious shifts in the mid-infrared to FIR colours 
with increasing galactocentric radius. 

Hi) The variation of the amplitude of the infrared brightnesses 
with grain size for a fixed galaxian position (radiation fields). 

As apparent from Fig. [7] the amplitude of the infrared brightness 
increases with increasing grain size, with a small increase for 
the small grains and a bigger increase for the bigger grains. 
Since the width of the SEDs is approximately constant (for 



18 Cristina C. Popescu et al.: Modelling the spectral energy distribution of galaxies. 




10 100 1000 10 100 1000 10 100 1000 



Fig. 7. Infrared brightnesses for grains of different sizes (plotted as different curves in each panel) and various composition: Si (left 
panels), Gra (middle panels) and PAH + (right panels), heated by the diffuse radiation fields calculated for the best fit model of our 
prototype galaxy NGC 891. Infrared brightnesses for PAH are not plotted in this figure. The colour coding is as follows: red is for 
grains with radius a < 0.001 fan (0.00035, 0.00040, 0.00050, 0.00063 and 0.00100 fan), green is for grains with 0.001 < a < 0.01 
(0.00158, 0.00251, 0.00398, 0.00631, 0.01000 fan) and blue is for grains with a > 0.01 (0.0316, 0.10000, 0.31623, 0.7943 /mi). The 
grain sizes considered in this plot are the same as those plotted in Fig. [6] Going from the top to the bottom panels the calculations are 
done for different positions in the model galaxy: R = Opc, z = Opc; R = 15000pc, z = Opc; R = Opc, z = 600 pc; and R = 15000 pc, 
z = 600pc. 



a constant set of parameters), the amplitude of the infrared 
brightness will scale with the area under their SEDs, namely 
with the energy absorbed (per unit a 2 ). For a given radiation 
field u m d, Eq. 1331 tells us that this is determined only by the 
optical properties of the grains in the UV and optical regime. 
Indeed, overall the Q a i, s increases with increasing grain size, 
with a smaller trend for smaller sizes, and a bigger trend for 
bigger grains. Especially in the case of silicate, the big grains 
show a tendency for higher efficiency in absorbing optical 
photons, which will boost the amplitude of their SEDs due to a 
higher proportion of red photons being absorbed. Thus, bigger 
grains will absorb more red photons than the smaller grains, 



which, for a fixed colour of the radiation fields, will allow big 
grains to have a higher increase in their infrared brightnesses 
with increasing grain size. 

iv) The variation of the amplitude of the infrared brightnesses 
with grain size for a variable galaxian position ( radiation fields). 

As mentioned in Hi), the bigger grains have a faster increase in 
the amplitude of their infrared SEDs with increasing size than 
the small grains, due to the increase in the efficiency of absorb- 
ing red photons. If the radiation fields will also change in colour, 
due to their spatial variation, this will induce an additional differ- 
ence in the amplitude of the SEDs of big and small grains. Thus, 
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Fig. 8. The best fit model SED of NGC891, corresponding to r£ = 3.5, SFR = 2.88 M /yr, F = 0.41, together with the observed data 
in the MIR/FIR/submm (plotted as rhombus symbols with error bars). The observed flux densities and the corresponding references 
are given in Table Q] The calculated total intrinsic stellar SED is plotted as black line and rhombus symbols. The predicted stellar 
SED given by population synthesis models is plotted as dark blue line. The intrinsic stellar SED of the young stellar population in 
the diffuse component (the fraction 1 — F escaping the star-forming complexes) is plotted with mauve line. The intrinsic SED (only 
diffuse) of the old stellar population, as derived from the optimisation of the optical/NIR images of NGC891, is plotted with light 
blue line. The corresponding SED for the bulge is plotted with pink line. The different symbols overplotted on the UV/optical SEDs 
indicate the wavelength at which the radiative transfer calculations were performed. The predicted PAH and dust emission SED is 
given with red line. Also plotted are the two main components of infrared emission: the predicted diffuse SED (dark-green line) and 
the predicted clumpy component associated with the star-forming regions (light-green line). 



at small galactocentric radii, where the radiation fields are red- 
der (see Fig. H]i, more optical photons are available to boost the 
amplitude of the SEDs of the big grains. The trend of increas- 
ing the relative contribution of big grain emission to the small 
grain emission with decreasing galactocentric radius is apparent 
from Fig. [7] especially for the case of silicates, which have a 
broad range in grain sizes. Obviously the increase in the inten- 
sity of the radiation fields with decreasing galactocentric radius 
will also increase the amplitude of the SEDs, but this will pro- 
duce an overall boost for both small and big grains. 

To conclude, the increase in the relative contribution of big 
grain emission to the small grain emission is a strong function 
of the colour of the radiation fields, and, unlike the wavelength 
dependence, does not depend on the heating mechanism of the 
grains. This also shows that models that assume a fixed colour 
of the radiation fields (e.g. that of the local interstellar radiation 
fields) will incur systematic errors in the mid-infrared to FIR 
colours with increasing galactocentric radius. 

Finally, the increase in the amplitude of the 8 fim PAH fea- 
tures with respect to the 3.3 yum feature at small galactocentric 
radius is also due to the redder radiation fields, which will pro- 
vide additional optical photons capable of exciting the vibra- 
tional transition at around 8 pm, but not energetic enough to ex- 
cite the transitions around 3.3 /urn. 

3.4. The integrated SED 

The best fit SED, obtained by spatially integrating the solution 
described above for the diffuse dust emission component and 



adding it to the solution for the dust emission from the localised 
component from the star-formation regions, is shown in Fig. [8] 
together with the observed FIR/submm data that were used to 
constrain the model solution. Details on the observed flux den- 
sities used in the plot are given in Table [2] The best fit solution 
corresponds to T f B = 3.5, SFR = 2.88M /yr and F = 0.41. 

Compared with the solution obtained in Paper I {t b = 4.1, 
SFR = 3.8M /yr and F = 0.22) there are changes in all three 
free parameters. The decrease in t b is solely attributed to the 
fact that we used improved observational data to fit the submm 
points. If we had used the original data the solution for opacity 
would be unchanged. The decrease in S FR and the increase in F 
is essentially due to the combination of adding the contribution 
of the young stellar population in the optical, and of including 
the wavelength dependence of the fraction of photons escaping 
the star-forming regions. 

With the red line we show the model fit for the total dust 
emission of NGC 891. We also show the main components 
of the dust emission, the diffuse component (dark-green line) 
and the localised emission from the clumpy component (light- 
green line). One can see that overall most of the dust emis- 
sion is powered by the diffuse component. Our solution gives 
for the total dust luminosity L'°'" s l t = 9.94 x 10 36 W, of which 

L dtJt = 689 x !0 36 W is emitted in the diffuse medium (69%). 
From the figure it is apparent that the diffuse component domi- 
nates the emission longwards of 60 /jm and shortwards of 20 fim. 
Thus, most of the emission in the FIR and in the NIR/MIR 
(PAH region) is diffuse. It is only at intermediate wavelengths 
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Table 2. The total MIR/FIR/submm flux densities of NGC 891. The observed flux densities available only within a fixed aperture 
were extrapolated to the whole galaxy using the best fit model of NGC 891. 



wavelength 




error(F,,) 


telescope 


reference 




/im 


Jy 


Jy 








12 


5.46 


0.82 


IRAS 


Rice et al. (1988); Sanders et al. (2003) 


average 


25 


7.39 


1.11 


IRAS 


Rice et al. (1988); Sanders et al. (2003) 


average 


60 


63.8 


9.6 


IRAS 


Rice et al. (1988); Sanders et al. (2003) 


average 


100 


185 


28 


IRAS 


Rice et al. (1988); Sanders et al. (2003) 


average 


170 


193 


18 


ISO 


Popescu et al. (2004) 




200 


165 


17 


ISO 


Popescu et al. (2004) 




450 


40.95 


8.0 


SCUBA 


Israel et al. (1999) 


extrapolated 


850 


4.85 


0.6 


SCUBA 


Israel et al. (1999) 


extrapolated 


1200 


0.92 


0.13 


IRAM 


Guelin et al. (1993) 


extrapolated 



(20 - 60 i-im) where the localised dust emission within the star- 
forming complexes dominates. In the IRAS 25 yum and Spitzer 
24 yimi bands almost all the emission is predicted to come from 
the star-forming complexes, suggesting the efficacy of this band 
as a direct tracer of SFR. Indeed empirically studies have shown 
that of the Spitzer bands it is the 24ji/m band which is more 
closely related to star-formation (e.g. Calzetti et al. 2010). 

In Fig. [8] we also show the corresponding intrinsic stellar 
SED (as would be seen in the absence of dust) of NGC 891, 
together with the stellar emissivity components. In all cases the 
symbols indicate the wavelengths at which the radiative trans- 
fer calculations were performed. The black line represents the 
total stellar emission produce in the galaxy. In the UV range 
this is given by the population synthesis models (plotted as the 
blue line), as explained in Sect. 12.3.31 A fraction of this emis- 
sion is locally absorbed in the star-forming complexes, while the 
remaining 1 — F escapes in the diffuse young stellar disk (the 
mauve line). Thus, the difference between the black line and the 
mauve line in the UV represents the localised (and wavelength 
dependent) absorption of stellar light (see Eq.[2Th. In the optical 
range the total emission is given as a sum of the emission from 
the young stellar disk (mauve line plots only the diffuse emis- 
sion from the young stellar disk, which is slightly different from 
the total emission of the young stellar disk, due to the local ab- 
sorption), the old stellar disk (light blue line) and the bulge (pink 
line). 

One important result of such calculations is to determine the 
fractional contribution of the different stellar components to the 
dust heating. For the case of NGC 891 we derive the following 
fractions: 11% for bulge, 20% for the old stellar disk, 38% for 
the young stellar disk and 31% for the star forming complexes. 
This means that the young stellar populations are responsible for 
69% of the dust heating while the old stellar populations account 
for the remaining 31%. 



4. Fidelity of the model 

In order to fit observed dust/PAH emission SEDs of real galaxies 
like NGC 891 it was necessary to make some basic assumptions 
and approximations, which, however, have implications not only 
for the SED of the integrated re-radiated light, but also for the ge- 
ometrical characteristics of this light and of the direct UV/optical 
stellar light. The main assumptions are the existence of a diffuse 
dust component associated with the young stellar population, the 
approximation of this dust in the form of an exponential disk and 
the utilisation of a fixed spectral emissivity law for the young 
and old stellar populations. Here we check whether our model 
is consistent with the available observational constraints beyond 



those of integrated dust emission SEDs and evaluate the limita- 
tions imposed by the approximations. 

4.1. The existence of a diffuse dust component associated 
with the young stellar population 

A fundamental aspect of our model was the inclusion of a sec- 
ond disk of dust associated with the young stellar population, 
which was taken to mimic the diffuse dust that is known to ex- 
ist in spiral arms through direct observations in the FIR, both in 
the Milky Way and in nearby well resolved galaxies. The second 
dust disk was originally introduced in Paper I to provide the ob- 
served level of submm emission in the spatially integrated SED. 
As discussed in Paper I and in Popescu & Tuffs (2005), the ad- 
ditional quantity of dust needed to fit the observed submm data 
cannot be provided by clumpy optically thick dust, whether such 
dust is in star formation regions or in passive quiescent clumps. 
This conclusion has been reinforced by the utilisation of the im- 
proved dust emission templates for star formation regions intro- 
duced in Sect. I2.8I since these are now empirically constrained 
by submm and near-mm data. 

More fundamentally, although one is in principle free to add 
dust to the model in any form one likes to reproduce the observed 
level of emission deep in the submm, if this dust is self-shielded, 
it will in practice struggle to supply the necessary luminosity to 
fit the FIR flux density peak of spiral galaxies at around 160 fim. 
To peak at around 160 micron the dust grains must be heated by 
radiation fields at around 1 Habing, which is indeed the illumina- 
tion of the diffuse dust. Heating at around 100 Habing (as at the 
PDR surface in embedded star-formation regions) would provide 
the luminosity but with too blue dust FIR colours, while the dust 
emission from the self-shielded dust would be too red as well 
as not providing the luminosity^- Of course, one solution would 
be to invoke star-forming complexes extending to radii of around 
10 times their actual sizes into the diffuse component, so that the 
dust illumination is reduced to the required levels. However this 
is akin to spreading the dust around the spiral arms, which is ex- 
actly what our second dust disk solution tries to mimic, so the 
difference to a diffuse dust component then becomes semantic. 

12 A possible way out of this conundrum might be to invoke heating of 
the self-shielded grains by the absorption of secondary NIR/MIR pho- 
tons emitted at the PDR surface for the star-formation regions. However, 
even if this process could provide enough luminosity to account for the 
160yum peak in the integrated SED of galaxies, it could not simultane- 
ously boost the MIR PAH emission, which require UV or blue optical 
photons for excitation. By contrast the integrated MIR PAH emission 
is known to be statistically related to the 160 pm emission from spiral 
galaxies (eg Bendo et al. 2008), typically accounting for ~ 15% of the 
total re-radiated starlight. 
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In general, the fit to data provided by our two-dust disk 
model is one of relatively large opacity for the central regions 
of local universe spiral galaxies (r g = 3.5-4; Paper I; Driver et 
al. 2007; see also Sect. 13.4b . with most of the opacity provided 
by the second dust disk. This is in contrast to the solution of low 
opacity (r^ = 1) and one disk of dust (obtained by Xilouris et 
al. 1997 and Bianchi 2007 from modelling the optical data only) 
which fails by a factor of ~ 3 to reproduce the submm emission 
(Popescu et al. 2000; see also Baes et al. 2010). 

Apart from the need to fit the submm data, and the physical 
considerations already outlined in Sect. 12.11 linking the second 
dust disk with the corresponding star-forming disk of gas, there 
is further empirical evidence for the existence of a second dust 
disk, as derived from comparison of the model predictions with 
other data, as follows: 

i) Comparison of model predictions for the stellar emissivity in 
the optical with population synthesis models 

We checked whether the predicted intrinsic SED of the old 
stellar populations in NGC 891, as derived from optimisation 
of the optical data only (Xilouris et al. 1999), together with 
the corresponding SFR needed to produce the total luminosity 
emitted by the dust, could be fitted by the population synthesis 
models of Kotulla et al. (2009). In Fig. [8] one can see that 
the predictions for the optical emission from the old stellar 
population (light blue line) in the B, V, and I bands fall severely 
below the predictions from the population synthesis model 
(dark blue line). In this plot we already included the second dust 
disk. In its absence the SFR needed to reproduce the energy 
emitted in the infrared would be even greater, so the discrepancy 
would increase. It is clear from here that there is a need for 
extra stellar luminosity presumably hidden by extra dust (the 
second dust-disk in our formulation). In Popescu et al. (2000) 
we assumed that the young stellar population was only emitting 
in the UV, while in the optical we only had the contribution from 
the old stellar populations, as derived by Xilouris et al. (1999). 
This assumption is not supported by the predictions of the 
population synthesis models, and indicates that the optimisation 
of the optical images can only reveal information about the old 
stellar populations and associated dust, but completely hides the 
information about the young stellar populations and associated 
dust. 

ii) Comparison of model prediction for the spatial distribution 
of PAH emission with observations 

Since our model with the second disk of dust and stellar emis- 
sivity was only fitted to the spatially integrated dust emission, 
it is possible to use the spatial information provided by the 
observations to check the predictions of the model. We have 
already used the predictions from Popescu et al. (2000) to 
compare the predictions of the model for the spatial distribution 
of infrared emission with the ISOPHOT data of NGC 891 at 
170 and 200 pm (Popescu et al. 2004). The comparison was 
done for the radial profiles, as the observed ISOPHOT images 
were only resolved in radial direction. The model predictions 
were found to be in excellent agreement with the observations. 
More recently Spitzer images of NGC 891 became available, 
which have a linear resolution of lOOpc at 8jt/m, comparable 
to the thickness of the second dust disk. Since we have now 
included the PAH features in the model, which dominate the 
8 pm emission, these data can be used to test the predictions for 
the vertical distribution of PAH emission in NGC 891. For this 
comparison we first subtracted the component of direct stellar 
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Fig. 9. Comparison between the 8 pm Spitzer vertical profile of 
NGC 891 integrated over longitude (dashed line) and our corre- 
sponding model predictions (solid line). Before comparison with 
the model, the 8 pm image had the stellar component subtracted 
off. 



light at 8 pm (estimated by multiplying a Spitzer 3.6 pm image 
by a factor of 0.26; Helou et al. 2004). Fig. [9] shows again an 
excellent agreement. Thus the model containing a second dust 
disk can predict both the overall level of emission and the shape 
of the predicted vertical profile. 

Hi) Comparison of model predictions for the attenuation- 
inclination relation with observations 

A strong independent evidence for the existence of a second 
dust disk came from the attenuation-inclination relation. It was 
shown in Popescu & Tuffs (2009) that while a two dust-disk 
model with higher central face-on opacity can reproduce the ob- 
served attenuation-inclination relation, a single dust disk model 
with t£ = 1.0 completely fails to reproduce the observed data. 
The attenuation-inclination relation is an especially sensitive 
test, as the rise in attenuation with inclination will very strongly 
depend on the relative scaleheights of the assumed dust layers 
and stellar populations. In particular it is an independent test for 
the existence of the second component of dust represented by 
the second dust disk. We should also mention here that one of 
the strength of this test is that it is completely independent of the 
assumed dust emission properties. Thus, this would also seem 
to rule out a one-dust disk model with low opacity but with dust 
grains having modified optical properties in the submm (e.g. 
enhanced submm emissivity), as proposed by Alton et al. (2004) 
and Dasyra et al. (2005). 

iv) Comparison of model predictions for the vertical profiles in 
the optical and NIR in NGC 891 with observations 
The solution with a second dust-disk and extra luminosity com- 
ing from the young stellar population emitting in the optical 
bands was further used to predict the appearance of the galaxy 
in the optical bands. In particular the average vertical profiles 
have been used to test whether the dip produced in the plane 
of the galaxy by the dust layer has the right deepness. Fig. [10] 
shows a comparison between the I band averaged vertical pro- 
file obtained from the observed images of NGC 891 (Xilouris 
et al. 1998) and the corresponding model predictions for three 
cases. The first case (left panel) is the original solution obtained 
by Xilouris et al. (1999), which only includes the old stellar disk 
and associated dust. The second case (middle panel) is our two 
dust-disk model, but without the inclusion of the young stellar 
disk. Finally, the right panel shows again the predictions for our 
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Fig. 10. A comparison between the observed averaged vertical profile of NGC 891 in the I band (solid line) and the corresponding 
model predictions (dashed line) for three cases: the one dust disk model (left panel) with r B — 1 .0 and one stellar disk (the old stellar 
disk); the two dust-disk model (middle panel) with = 3.5 and one stellar disk (the old stellar disk); the two dust-disk model (right 
panel) with t b — 3.5 and two stellar disks (the old and the young stellar disk). 



two-dust disk model, but with the inclusion of extra luminosity 
coming from the young stellar disk. One can immediately see 
that the inclusion of a second dust disk which is not accompa- 
nied by a stellar luminosity component would produce a stronger 
dust lane (a larger depth in the vertical profile) than shown in the 
observed images. This is a clear indication that there is a need 
for extra luminosity, as also indicated by population synthesis 
models (see i) above). In fact the two dust disk model with two 
(old and young) stellar components is able to reproduce the ob- 
served data better than the one disk model, especially in the B 
band. 

In the K band the model with the two dust disks and two 
stellar components predicts a somewhat more prominent dust 
lane than observed. This would indicate either a shorter scale- 
length for the second dust disk (though this would be difficult to 
reconcile with the excellent fit we found in Paper I for the ra- 
dial profile at 850 /mi), or more luminosity in the young stellar 
disk than predicted by the population synthesis model. It was al- 
ready noted by Dasyra et al. (2005), Bianchi (2007) and Bianchi 
(2008) that the second dust disk shows its maximum effect in 
the K band, since it is only there that the first disk of dust be- 
comes transparent and is therefore not shielding the second disk 
of dust. A possible alternative reason for the difficulty in fit- 
ting the vertical profile in K-band is that, due to its very small 
scale height, the appearance of the second dust disk at these 
long wavelengths would be easily blurred in real galaxies if per- 
turbations from co-planarity occur, even if any such perturba- 
tions had relatively small amplitudes. Although the scaleheight 
of the molecular layer in Milky Way is around 90 pc, as adopted 
in our model for the second dust layer, it is well known that 
CO and other tracers of star-formation exhibit systematic verti- 
cal displacements from the mean place known as "corrugations" 
(Spicker & Feitzinger 1986 and references therein). The discov- 
ery of Matthews & Uson (2008) of a non-planar disk in star- 
formation tracers in an isolated galaxy other than the Milky Way, 
is a strong evidence that non-planarity is a rather frequent phe- 
nomena. Matthews & Uson found that undulated patterns with 
amplitude of ~ 250 pc are visible in particular in the distribution 
of the young stellar population and the dust, suggesting that the 
process leading to the vertical displacements may be linked with 
the regulation of star formation in galaxies. The effect of corru- 
gations will mean that the edge-on thin disk will not be seen as 
such in the optical, but only the distribution corresponding to the 
amplitude of the corrugations. This would be more than enough 
to blur and hide any dust layer in the NIR. Nevertheless, the 



young stars will be still perfectly correlated with the thin layer 
of dust and gas with a (local) scaleheight of 90 pc, since the effect 
of corrugations is purely gravitational. So the basic vertical geo- 
metrical coupling between stars and dust which gives rise to the 
dust emission from the thin disk, as calculated in this model, still 
applies. To conclude, if corrugations occur, the second dust-disk 
will tend to be blurred, making dust lanes in edge-on galaxies in 
the NIR less prominent. 

4.2. Approximating the spiral arm component with an 
exponential disk 

Even if it is accepted that there is a diffuse dust layer associ- 
ated with the young stellar population, we must still evaluate the 
effect of artificially distributing this stellar emissivity and dust 
opacity in an exponential disk instead of a spiral arm pattern. It 
is clear that this will completely prohibit a comparison of surface 
brightness distributions in face-on systems, but here we are only 
concerned with the spatially integrated dust emission SEDs and 
the effects of this approximation on these. 

To do this we ran a simulation with the parameters corre- 
sponding to the best fit solution of NGC 891 (see Sect. 13.4b . i.e. 
we kept the same luminosity for the young stellar population and 
the same amount of dust, but we redistributed the correspond- 
ing stellar emissivities and dust opacities in a spiral pattern. We 
modeled the spiral pattern using 3 circular arms, with a radial 
distribution as given in Fig. QT| The vertical distribution is the 
same as for the case of the thin exponential disk. This means 
that for any given line of sight through the spiral arm the opac- 
ity is higher than for the case of a second dust disk, and for any 
given line of sight through the interarm regions the opacity is 
lower than for the case of the second dust disk - being just the 
opacity of the first dust disk. So we will have a solution with a 
high contrast between arm and interarm regions, that can be still 
characterised by the same central face-on opacity - which is the 
effective central face-on opacity if all the dust were distributed 
in an exponential disk instead of in a spiral arm. 

After calculating the radiation fields and the infrared emis- 
sivities using the same procedure outlined in Sect [2] and [3] we 
obtain a total integrated infrared SED that looks very similar to 
that obtained for the exponential disk case (see Fig.[T2l. The in- 
tegrated dust luminosity is only 5.5% lower than for the standard 
model. Here we should mention that the attenuation-inclination 
curve (not plotted in this paper) obtained for the solution with 
a spiral pattern is also almost identical to that obtained for the 
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Fig. 11. The radial distribution used to test the effect of a spi- 
ral pattern for the young stars and associated diffuse dust (solid 
line), overplotted on the distribution for the standard model with 
an exponential disk (dotted line). The radial coordinate is given 
in units of exponential scalelengths, where the scalelength of the 
exponential disk is h tAlsk = 5670 pc. 
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Fig. 12. Comparison between model solutions for the dust emis- 
sion SEDs calculated under the assumption that the emissivity 
of the young stellar population and the opacity of the associated 
dust are distributed in a thin disk (solid line) and in a spiral pat- 
tern (dashed line), respectively. The solution is for the best fit of 
NGC891. 



case of an exponential disk. The insensitivity of the shape of the 
attenuation-inclination curve to the inclusion of a spiral pattern 
was already demonstrated by Semionov et al. (2006). In passing, 
we also note that Misiriotis et al. (2000) showed that the appear- 
ance of simulated dusty spiral galaxies seen edge-on, calculated 
using a spiral structure, does not differ from that calculated using 
a pure exponential disk. 

The similarity in the solutions obtained for both the inte- 
grated infrared SEDs and for the attenuation-inclination relation 
reassures us that the approximation of a second exponential thin 
disk of stellar emissivity and dust opacity is an excellent one 
when making predictions for the spatially integrated SEDs. This 
is in contrast to the large change in the shape and amplitude of 
the infrared SEDs (as we will show in Sect. [5]l and attenuations 
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Fig. 13. Comparison between model solutions for dust emission 
SEDs calculated using intrinsic SEDs of the old stellar popu- 
lation peaking in the J (solid line) and K band (dashed line), 
respectively. The solution is for the best fit of NGC 891. 



of stellar light (see Fig. IB. II ) when changing the main free pa- 
rameters of our model. 

4.3. The approximation of a fixed spectral shape for the 
SEDs of the old and young stellar populations 

As discussed in Sect. 12.31 the dust emission SEDs are calculated 
by dividing the stellar emission SED into two fixed spectral tem- 
plates with differing spatial distributions - one UV-dominated 
corresponding to the young stellar population, weighted by the 
model parameter SFR, and one optically dominated, corre- 
sponding to the old population and weighted by the model pa- 
rameter old. This is a radically different approach to the handling 
of stellar emissivity SEDs compared to previous models for the 
panchromatic UV-submm emission of galaxies. It enables us to 
obtain the same solution for the dust emission, and therefore the 
reddening of any given galaxy, without the need to input trial 
population synthesis solutions for the full UV/optical/NIR stellar 
emission SED to the calculation of the dust emission. This new 
approach requires however that the solution to the dust emission 
is invarient to the assumed shape of the stellar emissivity SED 
within each of the two templates. Here we test this assumption 
in turn for the old and the young stellar emission templates. 

For the old stellar population the spectral shape was empir- 
ically derived from fitting the optical images of NGC 891. The 
spectral shape was consistent with the SED having a flux den- 
sity per unit frequency peaking in the J band. Reasonable varia- 
tions from this shape are to consider spectral templates peaking 
towards longer wavelengths, e.g. in the K band. We therefore 
consider a calculation in which the spectral shape was altered 
to allow for a brightening of the K band luminosity by 58% 
and a corresponding dimming of the J band luminosity by 48%, 
change that preserves the spectral integrated luminosity of the 
old stellar population. The result of this calculation from Fig.[T3l 
shows that the predicted dust emission SEDs are completely in- 
sensitive to (possible) changes in the spectral shape of the old 
stellar populations, as indicated from the overlap of the infrared 
SEDs. 

For the young stellar population the spectral shape was fixed 
from a combination of population synthesis models and empiri- 
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Fig. 14. Comparison between model solutions for dust emission 
SEDs calculated using different spectral templates for the intrin- 
sic SEDs of the young stellar population in the UV: our standard 
model (green line) and a model with bluer colours (red line), re- 
spectively. The solution is for the best fit of NGC 891. In both 
cases the the total luminosity of the young stellar population was 
kept fixed. 



cal fitting of the optical images. This resulted in a spectral tem- 
plate having a flux density per unit frequency peaking in the 
I band. Reasonable variations from this shape are to consider 
spectral templates peaking towards shorter wavelengths, e.g. in 
the V band. We therefore consider a calculation in which the 
spectral shape was altered to allow for a brightening of the V 
band luminosity by 34% and a dimming of the I band luminos- 
ity by 25%. As before this change was chosen to preserve the 
spectral integrated luminosity of the young stellar population. 
The resulting dust emission SED (not shown in this paper) is 
essentially indistinguishable from the one calculated using our 
standard model. There is only an increase in the predicted dust 
luminosity by 0.2%, which is to be expected due to the shift in 
the spectral peak of the young stellar population to shorter wave- 
lengths, which resulted in a larger fraction of stellar photons 
being absorbed by dust. So changing the colours of the young 
stellar population in the optical by as much as 59% produces 
changes of less of a percent for the predicted dust emission lu- 
minosity. 
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Fig. 15. Comparison between model solutions for diffuse dust 
emission SEDs calculated using different contributions of the 
young and old stellar populations, but preserving the total stel- 
lar luminosity. The solid line is the solution corresponding to 
the best fit of NGC891, but without the bulge contribution. The 
dashed line corresponds to an increase in the luminosity of the 
young stellar population in the diffuse stellar disk by a factor of 2 
and a decrease in the luminosity of the old disk stellar population 
by a factor of 1.68. 



Finally, it is important to check what the effect of changing 
the colours of the young stellar population in the UV is, since 
there it is these photons that make a significant contribution to 
the dust heating. We therefore made a more drastic change, by 
changing the slope of the UV SED from nearly flat (in F v ) to a 
monotonically decreasing slope between the 1350A and 2800 A 
spectral sampling points. We therefore brightened the 1350A 
data point by 24% and dimmed the 2800A flux by 38%, thus pro- 
ducing a very blue spectrum. As before we kept the overall lu- 
minosity of the young stellar population constant. The resulting 
dust emission SED is shown in the two panels of Fig. [14] plotted 
in a linear scale and for different cuts in luminosity to allow a 
better visualisation of the small differences in the SEDs. Again, 
since a bluer stellar SED was considered, this resulted in more 
stellar photons being absorbed by dust and a larger predicted 
dust luminosity, as seen from the plots. The overall increase in 
the dust luminosity was 1.7%, which is still a minor variation 
taking into account the dramatic change by 62% in the colours 
of the UV stellar SED. In fact this is almost surprising, but we 
should keep in mind that the bluer UV photons have higher prob- 
ability of absorption in the star forming complexes, and therefore 
the escaping radiation illuminating the diffuse component would 
still be strongly reddened due to the local absorption. 

By contrast we also did a calculation in which we changed 
the relative contribution of the young to old stellar populations 
and kept fixed the total stellar luminosity. We have used the so- 
lution for NGC 891 for the diffuse component, set the bulge- 
to-disk ratio to and used this as a reference SED. We then in- 
creased the luminosity of the young stellar population by a factor 
of 2 and decreased the luminosity of the old stellar populations 
in the disk by a factor of 1.68, which preserves the total stel- 
lar luminosity. The resulting SEDs are shown in Fig. [15] This 
time there is a significant change in the SEDs, followed by an 
increase in the total dust luminosity by 33% and an increase in 
the MIR emission by a factor of ~ 2. This shows again that the 
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main factors influencing the dust emission SEDs are the overall 
luminosities of the old and young stellar populations, and not 
(possible) variations in the spectral shape of the template stellar 
SEDs. 

4.4. The relative opacity of the first and second dust disk 

A further approximation of the model is the fixed ratio between 
the opacity of the first and second dust disks. While most of 
the fixed parameters have been calibrated to empirical relations, 
we acknowledge that this ratio has been only calibrated to the 
value of NGC 891 (close to NGC 5907). The validity of this as- 
sumption has been succesfully tested in a statistical sense on the 
attenuation-inclination relation in Driver et al. (2007). However 
this assumption still needs to be proven on an object-by-object 
basis by fitting the panchromatic SEDs of galaxies. As we will 
show in Appendix [D] we have identified a parametric test to 
potentially flag out galaxies that may not follow the colour- 
luminosity relation predicted by our model due to other geomet- 
rical characteristics. 

5. The library of SEDs for the diffuse component 

We have created a library of diffuse SEDs that spans the param- 
eter space of 4 parameters: tL, S FR', old, and With the 
exception of S FR', these are the main parameters of the model. 
S FR' is defined in terms of the primary parameters S FR and F: 

SFR' = SFR x( l - F) (45) 

which is equivalent with the S FR for the case that F — 0. This 
definition is possible because in our model we made the assump- 
tion that the wavelength dependence of the fraction of escape 
photons from the clumpy component into the diffuse one is fixed, 
allowing us to separately calculate the SEDs for the diffuse com- 
ponent and for the clumpy component. In this formulation S FR' 
is an effective star formation rate powering the diffuse dust emis- 
sion. 

In Table [3] we give the parameter values at which the model 
for the diffuse component has been sampled. In total we have 7 
(for t{) x 9 (for SFR') x 9 (for old) x 5 (for B/D) = 2835 com- 
binations, corresponding to 2835 simulated SEDs. All the model 
SEDs are available in electronic form. The choice of the param- 
eter values was done such that it covers the parameter space of 
local universe galaxies but also the asymptotic values of these 
parameters. Thus for r^j we considered the range [0.1,8.], which 
means galaxies with almost no dust at all to galaxies having over 
twice as much opacity as the average value found from the statis- 
tical analysis of the Millennium Galaxy Catalogue (Driver et al. 
2007). For S FR' we considered the range [0.,20.] M /yr, which 
means galaxies with no recent star-formation activity (no heating 
from the young stellar population) to galaxies having very high 
star-formation rate (typical of starburst galaxies which should 
be outside the range of spiral galaxies addressed by this paper). 
For the old stellar population the parameter old was scanned 
analogous to SFR. Finally, we considered galaxies spanning the 
whole Hubble sequence, from bulgeless galaxies {B/D = 0.0) to 
bulge dominated galaxies (B/D = 2.0). This choice of param- 
eters allows a smooth overlap with the parameter space of star- 
burst galaxies, very quiescent galaxies, early type galaxies, and 
higher redshift galaxies. It also allows rare (unexpected) cases 

13 We note that the fixed (and calibrated) parameters of the model are 
given in the Appendix [Ej 



Table 3. The parameter space sampled by our library of dust and 
PAH emission SEDs (for the diffuse component only). 





n 


SFR' 

M /yr 


old 


B/D 


1 


0.1 


0.0 


0.0 


0.00 


2 


0.3 


0.1 


0.1 


0.25 


3 


0.5 


0.2 


0.2 


0.50 


4 


1.0 


0.5 


0.5 


1.00 


5 


2.0 


1.0 


1.0 


2.00 


6 


4.0 


2.0 


2.0 




7 


8.0 


5.0 


5.0 




8 




10.0 


10.0 




9 




20.0 


20.0 





to be considered. As we will show in Sect. 16.51 the locus of the 
FIR colours corresponding to the parameter space defined by the 
values in Table [3] overlaps quite well with the locus of observed 
FIR colours of real life spiral galaxies. 

To compute the library of diffuse SEDs we first calculated 
a library of radiation fields, computed for the main stellar com- 
ponents of our model: young stellar disk, old stellar disk and 
bulge, for the 7 values of opacity used in Tableland for the 16 
UV-optical wavelengths detailed in Sect. 12.31 For the bulge and 
the old stellar disk only 6 optical wavelengths were used, as in 
our model we assume that these stellar components have neg- 
ligible emission in the UV range. In total we created a library 
of 196 data cubes of radiation fields, sampled at 22 radial po- 
sitions and 12 vertical positions (264 spatial points within the 
model galaxy). We then created the library of temperature dis- 
tributions, for the 4 grain compositions used in our model (sili- 
cate, graphite, PAH , PAH + ), the 7 values of opacity, 9 values of 
S FR', 9 values of old and 5 values of B/D. In total we computed 
a library of 1 1340 data cubes of temperature distributions, each 
sampled at 264 spatial points within the model galaxy and for 
all grain sizes contained in the dust model. For each of the tem- 
perature distribution data cubes we created corresponding cubes 
of infrared emissivity. In total a library of 1 1340 files of infrared 
emissivities were calculated. 

When used to fit observed panchromatic SEDs, the library 
of simulated dust and PAH emission SEDs should be used in 
conjunction with the library of simulated attenuations of stellar 
light recalculated in this paper, taking into account the B/D ratio 
of the galaxy and its inclination, as described in Appendix ID1 

6. Predicted variation of the dust emission SEDs 
with the main parameters of the model. 

6.1. Variation of the SEDs with r£ 

i) Amplitude of the SED 

As expected, the amplitude of the diffuse dust and PAH 
emission SEDs increases with increasing optical depth for the 
optical thin cases and tends to a saturation value for the optically 
thick cases. This is seen in Fig. [16] from the bunching of the 
blue and black curves on one hand (t b = 4, 8) and from the 
big gap between the green and red curves on the other hand 
(t b = 0.1,1). An interesting feature of the SEDs is the fact 
that the ratio between the FIR and MIR (PAH) amplitudes 
increases with increasing opacity for the models where the 
stellar luminosity has a higher contribution from the old stellar 
population with respect to the young stellar population (but 
where the young stellar population has still a non-negligible 
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Fig. 16. Integrated dust and PAH emission SEDs for the diffuse component, for model galaxies with different face-on opacities 
(plotted as different curves in each panel). All models shown in this figure are for pure disk galaxies (BID = 0.0). From left to 
right the panels show model galaxies with various levels of contribution from the old stellar populations (old = 0.1 (0.0), old = 2.0 
and old = 20.0). From top to bottom the panels show models with various levels of SFR (SFR' = 0.0, 2.0, 20.0 M /yr). The colour 
coding is as follows: green is for = 0.1, red is for r£ = 1 .0, blue is for = 4.0 and black is for = 8.0. 



contribution). This can be seen for example on the right bottom 
panel of Fig. [16] where the green curve (r^ = 0.1) shows 
almost identical levels in the FIR peak and in the PAH features, 
while the black curve (r^ = 8) shows two order of magnitude 
difference in the FIR to MIR levels. This change in the FIR 
to MIR colour is due to the fact that with increasing optical 
depth the disk becomes first optically thick to the UV radiation 
(provided by the young stellar population), while still being 
relatively transparent to the optical photons (mainly provided 
by the old stellar population). This will have the consequence 
that the PAH emission (and small grain emission), which is 
mainly heated by the UV photons will tend to a saturation level, 
as expected for the optically thick case, while the FIR emission 
will be still boosted by the optical photons, which have not 
reached the optically thick limit yet. Apart from this effect, 
which is simply related to the optical properties of the grains, 
there is an additional geometrical effect which will boost the 
FIR-to-MIR colour with increasing opacity. This is due to the 
fact that the old stellar populations have larger scale-heights 
than the young stellar populations, which will mean that a larger 
proportion of the optical photons will be less confined to the 
regions of higher optical depth, and will therefore provide extra 



heating in the optically thin regions of the galaxy. 



ii) Peak of the SED 



Another interesting feature of the plots is the fact that the 
peak of the SEDs from diffuse dust shifts towards longer wave- 
lengths with increasing opacity for the optically thick regime, 
for the models where the stellar luminosity has a higher contri- 
bution from the young stellar populations with respect to the old 
stellar population. This can be seen for example in the bottom 
row of Fig. [16] where the trend of shifting the peak of the SEDs 
becomes increasingly less strong in moving from the left to the 
right panel, where the solution changes from a young to an old 
stellar population dominance of the stellar SED. As before, when 
moving to a solution where more optical photons are available, 
these will provide extra heating to the dust, due to their more op- 
tically thin regime, both because of their optical properties and 
because of their spatial distribution (as discussed before), which 
in turn will cancel the shift towards cooler SEDs. 
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Fig. 17. Integrated dust and PAH emission SEDs for the diffuse component, for model galaxies with different levels of SFR (plotted 
as different curves in each panel). All models shown in this figure are for pure disk galaxies (B/D = 0.0). From left to right the panels 
show model galaxies with various levels of contribution from the old stellar populations (old = 0.1, old = 2.0 and old = 20.0). 
From top to bottom the panels show models with various face-on B band opacities (t^ = 0.1, = 1.0, and = 8.0). The colour 
coding is as follows: green is for SFR' = 0.0M o /yr, red is for SFR' = 2.OM /yr, blue is for SFR' = 10.0M o /yr and black is for 
SFR' = 20.0M Q /yr. 



6.2. Variation of the SEDs with SFR (young stellar 
population) 



The increase in the contribution of the young stellar population 
to the stellar SEDs will have the effect of increasing the MIR 
to FIR level of the diffuse dust emission SEDs, with essentially 
no change in the submm level (see Fig. fTTb . At the same time 
the peak of the dust emission SEDs will shift towards shorter 
wavelengths. Overall this will result in warmer SEDs, with both 
FIR SED peaks and MIR-to-FIR colours becoming systemati- 
cally bluer. It is interesting to note that the effect of increasing 
SFR is completely different from the effect of increasing r^, 
showing that the two parameters are completely orthogonal. One 
should also notice that the trend in bluer SEDs with increasing 
S FR becomes less pronounced for models with higher contri- 
bution from the old stellar populations to the stellar SEDs (see 
trends in moving from the left column to the right column of 

Fig. mi). 



6.3. Variation of the SEDs with old (the old stellar population) 

The increase in the luminosity of the old stellar population pro- 
duces a shift of the peak of the infrared SED of the diffuse 
dust towards shorter wavelengths (see Fig. [T8l ). similar to the 
case of increasing the luminosity of the young stellar population. 
However, the shift is accompanied by a decrease in the MIR-to- 
FIR ratio, opposite to the effect obtained in the case of increas- 
ing the luminosity of the young stellar population. Indeed, the 
MIR-to-FIR colours become cooler, due to the fact that the ad- 
ditional optical photons will mainly boost the FIR emission and 
not the PAH emission. Obviously this effect is less pronounced 
for models with higher luminosities of the young stellar pop- 
ulations (see trends from going from the left to the right hand 
column in Fig. [T8l and is enhanced for models with higher dust 
opacity (see trends for going from the top to the bottom row in 
Fig.[H. 

6.4. Variation of the SEDs with B/D (bulge-to-disk ratio) 

The variation in the bulge-to-disk ratio produces variations in the 
infrared SEDs with less dynamical range (see Fig. [T9b . simply 
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Fig. 18. Integrated dust and PAH emission SEDs for the diffuse component, for model galaxies with different contribution from the 
old stellar population (plotted as different curves in each panel). All models shown in this figure are for pure disk galaxies (B/D = 
0.0). From left to right the panels show model galaxies with various levels of SFR (SFR' = 0.0, SFR' = 2.0 and SFR' = 20.0). From 
top to bottom the panels show models with various face-on B band opacities (t^ = 0.1, tL = 1.0, and = 8.0). The colour coding 
is as follows: green is for old = 0.0, red is for old = 2.0, blue is for old = 10.0 and black is for old = 20.0. 



due to the smaller dynamical range in the values of the parame- 
ter B/D. Nevertheless, the peak of the SEDs is shifted to shorter 
infrared wavelengths with increasing B/D and the MIR-to-FIR 
colours get cooler, following the trends expected for a variation 
in the luminosity of an old stellar population. Thus the variation 
induced by the B/D ratio follows similar trends with the varia- 
tion induced by the old parameter. It is however expected that 
in most cases the value of the B/D to be inputted from observa- 
tions, and thus not to be a free parameter of the model. 

Overall it is interesting to observe that the 3 parameters of 
the model: r B , SFR and old are orthogonal parameters, produc- 
ing quite different effects in shaping the dust emission SEDs. 
For example it is clear that the only way to increase the level of 
the submm emission is through a variation of r£, as neither an in- 
crease in the luminosity of the old or of the young stellar popula- 
tions could account for this. In other words in the submm we are 
tracing dust column densities, while in the MIR and FIR we are 
tracing both dust column densities and heating sources. It is also 
obvious that if we want to have warmer MIR-to-FIR colours we 
have to increase the luminosity of the young stellar populations, 
while if we want to have cooler MIR-to-FIR colours we need to 
increase the luminosity of the old stellar populations, with differ- 



ent degrees of modulations due to variations in r£. Here we need 
to remember that for fitting total SEDs the template for the star- 
forming complexes must be added to the diffuse SEDs, boosting 
the MIR emission (see Sect. 16.51 below). 

6.5. Colour variation of the combined diffuse and localised 
dust emission SEDs 

Adding the emission from the localised dust emission template 
of star-formation complexes to the solution for the diffuse emis- 
sion to obtain the total emission adds extra variations to the 
dust emission SEDs, especially boosting the warm dust emis- 
sion from larger grains. This effect, controlled by the F parame- 
ter, is in practice rather orthogonal to the variations in the SED 
of the diffuse component discussed above, due to the different 
behaviour of the UV/optical attenuation as a function of F (see 
Appendix ICt and due to the fact that, over the parameter range 
of the model, the localised dust emission template is almost al- 
ways much warmer in the FIR than the predicted FIR colours 
of the diffuse dust emission. In fact the FIR colours of the dif- 
fuse dust emission only approach those of the PDR template for 
the extreme case S FR = 20 and old = 20 in an optically thin 
galaxy in which the UV illuminates the diffuse dust, as can be 
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Fig. 19. Integrated dust and PAH emission SEDs for the diffuse component, for model galaxies with different bulge-to-disk ratios 
(plotted as different curves in each panel). All models shown in this figure are for galaxies with a fixed contribution from the old 
stellar population (old = 2.0). From left to right the panels show model galaxies with various levels of SFR (SFR' = 0.0, SFR' = 2.0 
and SFR' = 20.0). From top to bottom the panels show models with various face-on B band opacities (r^ = 0.1, = 1.0, and 

= 8.0). The colour coding is as follows: green is for B/D = 0.0, red is for B/D = 0.5, blue is for B/D = 1.0 and black is for 
B/D = 2.0. 



seen by comparing the top RH panel of Fig. 18 with the template 
SED in Fig. 5. A further reason for the orthogonality of the tem- 
plate SED to the SED of the diffuse emission is that dust in the 
diffuse ISM is more efficient in chanelling absorbed UV energy 
into the PAH features than dust in the star-formation regions, 
due to the photo-destruction of PAH in star-formation regions 
(see Sect. 12.8b . Overall, solutions for the total infrared emission 
require most (though not all) of the emission in the 20 - 60 /mi 
range to be localised emission from star-formation regions, with 
the MIR and FIR/submm emission flanking this range being pre- 
dominantly diffuse in origin. 

The combined effect of the 5 main parameters of the model 
t^j, SFR, old, B/D and F will produce a large variation in 
the colours of the simulated SEDs. In Fig. [20] we plot the 
1 70/1 00 /mi colour versus the 100/60 /mi colour for our model 
SEDs. The SEDs were derived by combining all the simulated 
SEDs of the diffuse component from our library with the tem- 
plate SED for star-forming regions, whereas the combination 
was done for different values of the F factor (including the 
asymptotic values). This defines the locus in the colour-colour 
diagram occupied by our model SEDs. To check that the pa- 



rameter space covered by the models overlaps with the param- 
eter space of real life spiral galaxies, we overplotted the corre- 
sponding colours of the spiral galaxies with Hubble type earlier 
than Sd from the ISOPHOT Virgo Cluster Deep Survey (IVCDS; 
Tuffs et al. 2002a, 2002b; Popescu et al. 2002) and from the cat- 
alog of compact sources of the ISOPHOT Serendipity Survey 
(Stickel et al. 2000). One can see that except from a few outliers, 
the observed galaxies lie on the locus defined by our models. 
One can also observe that the position of NGC 891 in this dia- 
gram does not have any preferential place, but lies somewhere 
towards the very quiescent region of the models (as expected). 
More detailed comparisons with data and applications of the 
model will be given in a series of future papers. 

From Fig. [20] it is also obvious that our models embrace a 
somewhat larger area in parameter space than that defined by real 
life spiral galaxies. This was done on purpose, from the choice 
of the minimum and maximum values of the main parameters 
of the model. We tried to push the limited values towards the 
asymptotic values, in order to cover even rare (unexpected) cases 
that could occur in real life (see Sect. [5]and Table[3]for the range 
of parameter space). 
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Fig. 20. The 170/100 versus 100/60 micron colour-colour plot 
for our simulated SEDs (black diamonds). Overplotted with 
red triangles are the observed colours of the spiral galaxies 
from the ISOPHOT Virgo Cluster Deep Survey (IVCDS; Tuffs 
et al. 2002a, 2002b; Popescu et al. 2002) and the ISOPHOT 
Serendipity Survey (Stickel et al. 2000). Only spiral galaxies 
with Hubble type earlier than Sd are plotted. The white star rep- 
resents the position of NGC 891 in the colour-colour diagram. 

Most fundamentally, however, we emphasise that in the 
context of a RT solution it is the combination of the ob- 
served amplitude and colour of the dust emission SED with 
the self-consistently applied constraint on the attenuation in the 
UV/optical that leads to an inherently non-degenerate solution. 

7. Discussion 

In this paper we have assembled the components from which 
we can combine the predictions for the attenuation of starlight 
(as originally formulated in Paper III and re-expressed here in 
Appendix [B] and Appendix [Cb with the predictions for the SED 
of the dust/PAH re-radiated starlight (self-consistently calcu- 
lated as described in Sect. [2]) to extract the physical parameters of 
a star-forming galaxy from an observed UV/optical/NIR - mid- 
IR/far-IR/submm SED. 

The detailed mathematical prescription of how to combine 
the library of dust emission SEDs with the library of atten- 
uation to fit UV/optical to infrared/submm SEDs is given in 
Appendix iDl Here we discuss some of the issues regarding the 
applicability and physical relevance of our model. 

Firstly, our model is designed for data sets that have good 
enough optical images to derive B/D and i parameters. For the 
i parameter this should be the case for most optical surveys of 
local universe galaxies. For example, the SDSS, with an angu- 
lar resolution of 1.6", was able to derive inclination angles for 
galaxies with recession velocities ranging to * 10000 km/s. With 
modern surveys routinely having ss 0.5" resolution, inclinations 
should be available up to z ~ 0.1, which encompasses the local 
Universe galaxies. In cases where inclination angles cannot be 
derived, an expectation value may be assumed instead, and the 
quartile range of uncertainty in derived physical parameters re- 



sulting from a lack of knowledge of i can be derived by running 
separate optimisations for the quartile values of i. It should be 
noted that the inclination angle should not be considered as a 
free parameter, as its effect on the SED would not be orthogonal 
to the T f B parameter. The other parameter B/D should also be an 
input parameter for the routine and not a free parameter. In cases 
where information on B/D is missing, users are advised to use 
the correlation between B/D and apparent optical-NIR colours 
to estimate the B/D. 

The use in the model of the third parameter accessible 
through optical observations, 6 ga i, is qualitatively different to that 
of B/D and i. Whereas B/D and i can (for practical applications) 
only be used as constraints, 6 ga i may in practice be used either as 
a constraint (as in the example given above) or as a free param- 
eter. When Ogai is used as a free parameter, the angular surface 
area 8 2 gal effectively serves as free amplitude scaling factor for 
the amplitude of the diffuse dust emission (see Eq. ID.5l ). 

Secondly, our model is in its present form only designed to 
analyse star-forming galaxies. As noted above, when perform- 
ing an optimisation with the optical structural parameters fixed 
by direct observation, the optimisation no longer has any degree 
of freedom in terms of scaling parameters. This makes the model 
sensitive to a test of the fundamental assumption, that the dust is 
exclusively heated by stellar photons. The fits will be systemati- 
cally biased if dust is heated by photons from non-stellar sources 
like AGN or by some completely different channel. For example, 
by incorporating the optical and UV constraints, the model best 
fit parameters for SFR and t b will be likely to be skewed up- 
wards if an AGN heating dust with a harder UV photon spectrum 
than expected from a pure stellar source is present. In this case 
the dereddened UV/optical spectrum from step vii will however 
provide a direct flag for the presence of a bluer spectral emis- 
sivity in the UV than could be provided by stars. Potentially, 
therefore, our SED fitting technique can provide a method for 
recognising the presence of dust-obscured AGN activity. 

Thirdly, we note that the spectral form of the diffuse com- 
ponent of the dust emission SED given by Eq. ID.5l predicts that 
galaxies will lie on specific locii in a colour - surface brightness 
diagram for the dust emission, dependent on the surface density 
of UV/optical emissivity and t b . In cases where Q ga \ is known 
from optical measurements, any deviation of the observed posi- 
tions of galaxies in colour - surface brightness space from the 
positions predicted by the model can be therefore used as a non- 
parametric test of the fidelity of the geometry used in the model 
without having very specific a priori knowledge of the physi- 
cal properties of the observed sources. This would also serve as 
the most direct verification of the major role played by emission 
from diffusely distributed dust in our model. 

Finally, we caution the potential user of this model that, al- 
though RT solutions applied in this way are in principle a very 
powerful way to constrain physical properties of galaxies, their 
predictive power is ultimately reliant on an priori knowledge of 
the geometry of the system. Although we have specified this ge- 
ometry using empirical constraints derived from spiral galaxies 
of intermediate morphological type, and have endeavoured to ex- 
tend the applicability to galaxies of earlier or later types through 
the B/D parameter, constrained from optical imaging, the prac- 
tical range of morphological types where the advantages of the 
RT treatment outweigh deviations of real from assumed geom- 
etry still remain to be determined. This question can only be 
answered by statistical analysis, such as potentially realisable in 
a non-parametric way vis the colour-surface-brightness relation 
for the diffuse dust emission component (as proposed above) for 
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different morphological classes of galaxies, and will be the sub- 
ject of future studies. 

8. Summary 

In this paper we presented a radiative transfer model for spi- 
ral galaxies that self-consistently accounts for the attenuation of 
stellar light in the UV/optical and for the dust re-emission in the 
MIR/FIR/submm. We used this model to create a comprehen- 
sive library of dust and PAH emission SEDs and corresponding 
attenuations of stellar light (originally presented and described 
in Paper III and given in revised form in this paper) that can be 
used to routinely fit the panchromatic SEDs of large statistical 
samples of spiral galaxies. 

Our model has been calibrated for local Universe galaxies, 
and therefore its applicability should mainly lie within low red- 
shift galaxies. The models are also targeted to spiral galaxies. We 
have not attempted to model elliptical galaxies, starburst galax- 
ies, dwarf galaxies or AGN nuclei. 

The free parameters of our model for the calculation of the 
infrared SEDs are: 

- the central face-on opacity in the B-band r£, 

- the dumpiness factor F which defines the fraction of UV 
photons locally absorbed in star-formation regions and the 
corresponding fraction available to illuminate the diffuse 
dust, 

- the star-formation rate S FR, 

- the normalised luminosity of the old stellar population old, 

- the bulge-to-disk ratio B/D. 

For the parameter t b we have provided a relation for converting 
to total dust masses (see Eq.l44l>. 

The free parameters r£, F, and B/D are common to the set of 
parameters needed to predict the attenuation of the UV/optical 
stellar light (the fourth free parameter that affects the attenu- 
ation is the inclination i of the disk). Therefore the simulated 
dust emission SEDs presented here can be used in conjunction 
with the predictions for attenuation given in this paper for a self- 
consistent modelling of the UV/optical-IR/submm SEDs. Our 
model has also a parameter the angular size ga i, which can be 
used either as a constraint (fixed via optical images) or as a free 
parameter. When 6 ga i is used as a free parameter it determines 
the amplitude scaling of the fit. 

Since in our model we made the assumption that the wave- 
length dependence of the fraction of escape photons from the 
clumpy component into the diffuse one is fixed, we can sep- 
arately calculate the SEDs for the diffuse component and for 
the clumpy component. For the clumpy component we adopted 
the model of Groves et al. (2008), which we tuned to fit the to- 
tal MIR/FIR/submm SED of the ensemble of star-forming com- 
plexes in the Milky Way (Sect. 12.8b . For the diffuse component 
we defined an effective star-formation rate S FR' (corresponding 
to the total illumination by the young stellar population of the 
diffuse dust) and created a library of diffuse SEDs that spans the 
parameter space of tC, SFR', old, and B/D. In total we sam- 
pled 7 values in r^, 9 in SFR', 9 in old and 5 in B/D, making 
a total of 2835 simulated SEDs. This corresponds to a library of 
196 data cubes of radiation fields, sampled at 22 radial positions 
and 12 vertical positions (264 spatial points within the model 
galaxy), a library of 11340 data cubes of temperature distribu- 
tions and a library of 11340 files of infrared emissivities. We 
described (Appendix [Dj how in practice one can combine the 



predictions for the attenuation of starlight with the predictions 
for the SED of the dust/PAH re-radiated starlight to extract the 
physical parameters of a star-forming galaxy from an observed 
UV/optical/NIR - mid-IR/far-IR/submm SED. 

The analysis of the library of diffuse integrated SEDs 
(Sect. [6]) showed that the main parameters of the model are quite 
orthogonal, producing different effects in shaping the dust emis- 
sion SEDs. Specifically we have shown that: 

- The amplitude of the dust and PAH emission SEDs increases 
with increasing optical depth for the optical thin cases and 
tends to a saturation value for the optically thick cases. 

- The ratio between the FIR and MIR (PAH) amplitudes in- 
creases with increasing opacity for the models where the 
stellar luminosity has a higher contribution from the old stel- 
lar population with respect to the young stellar population. 

- The peak of the SEDs shifts towards longer wavelengths with 
increasing opacity for the optical thick regime, for the mod- 
els where the stellar luminosity has a higher contribution 
from the young stellar populations with respect to the old 
stellar population. 

- The increase in the contribution of the young stellar popu- 
lation to the stellar SEDs will have the effect of producing 
warmer SEDs, with both FIR SED peaks and MIR-to-FIR 
colours becoming systematically bluer. 

- The increase in the luminosity of the old stellar population 
produces SEDs with warmer FIR SED peaks but with redder 
MIR-to-FIR colours. 

- The variation in the bulge-to-disk ratio produces variations 
in the infrared SEDs with less dynamical range, simply due 
to the smaller dynamical range in the values of the parame- 
ter B/D. Overall the peak of the SEDs is shifted to shorter 
infrared wavelengths with increasing B/D and the MIR-to- 
FIR colours get cooler, following the trends expected for a 
variation in the luminosity of an old stellar population. 

Overall we found that the effect of increasing S FR is com- 
pletely different from the effect of increasing and works in the 
opposite direction from that produced by an increase in old. We 
also verified that the only way to markedly increase the level of 
the submm emission is through a variation of t b . 

The spatially integrated SEDs were derived from simulated 
images of dust emission, themselves calculated from 3D data 
cubes of infrared emissivity representing the response of grains 
(integrated over size distribution and composition) to the radia- 
tion fields derived from radiative transfer calculations. Our cal- 
culations take full account of the large variations in colour of the 
radiation fields as a function of position in the galaxy, which we 
show to be important (Sect. 13. U . Such variations lead to large 
trends in the grain temperature distributions (Sect. 13.2b used in 
the calculation of stochastic emission, leading to corresponding 
trends in the infrared brightnesses (Sect. [331 ) with position in the 
galaxy. In particular we concluded that: 

i) The shift of the peak of the infrared brightness as a func- 
tion of grain size strongly depends on the temperature of the 
grains, and therefore on the stochastic or non-stochastic nature 
of the heating mechanism. Since the heating mechanism depends 
both on grain size and on the intensity and colour of the radia- 
tion fields, it is clear that the shift cannot be described in terms 
of grain size only. This also shows that models that have a fixed 
grain size for the transition between the main heating mecha- 
nisms of dust, irrespective of the radiation fields, will lead to sys- 
tematic spurious shifts in the mid-infrared to FIR colours with 
increasing galactocentric radius. 
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ii) The increase in the relative contribution of big grain emis- 
sion to the small grain emission is a strong function of the colour 
of the radiation fields, and, unlike the wavelength dependence, 
does not depend on the heating mechanism of the grains. This 
also shows that models that assume a fixed colour of the radia- 
tion fields (e.g. that of the local interstellar radiation fields) will 
incur systematic errors in the mid-infrared to FIR colours with 
increasing galactocentric radius. 

Basis for the radiative transfer calculations were empirical 
constraints for the geometry of the large scale distributions of 
stellar emissivity and dust opacity of the translucent and opaque 
components of galaxies (Sect. l2~TT ). The translucent components 
were constrained using the results from the radiation transfer 
analysis of Xilouris et al., while the geometry of the optically 
thick components is constrained from physical considerations 
with a posteriori checks of the model predictions with obser- 
vational data. The main checks concern the nature of the diffuse 
distribution of dust associated with the young stellar population 
which is needed to predict the observed FIR/submm emission of 
real life galaxies. These checks (Sect.|U) showed that: 

i) Comparison of model predictions for the stellar emissivity 
in the optical with population synthesis models is not consistent 
with one dust disk model with one stellar disk component, and 
indicates the need for extra luminosity hidden by extra dust - the 
second dust disk and young stellar disk in our model formula- 
tion. 

ii) Comparison of model predictions for the vertical profiles 
of PAH emission with observations at lOOpc linear resolution 
reveals a good agreement, consistent with the existence of a thin 
dust layer, as represented by the second dust disk. 

iii) Comparison of model predictions for the attenuation- 
inclination relation with observations is also consistent with the 
existence of a thin dust layer, as represented by the second dust 
disk. 

iv) Comparison of model predictions for the optical surface 
brightness distributions of the edge-on galaxy NGC 891 with 
observations indicates that a two dust disk model with two (old 
and young) stellar components is able to reproduce the observed 
data better than the one disk model with one stellar disk. We also 
showed that a two dust disk model which is not accompanied 
by an extra stellar luminosity component produces a stronger 
dust lane than observed in real images. In the K band the model 
with the two dust disks and two stellar components predicts a 
somewhat more prominent dust lane than observed, indicating 
either a shorter scalelength for the dust disk (though this would 
be difficult to reconcile with the excellent fit we found in Paper I 
for the radial profile at 850 /mi) or more luminosity in the young 
stellar disk than predicted by the population synthesis models. A 
possible alternative reason for the difficulty in fitting the vertical 
profile in K-band is that, due to its very small scaleheight, the 
appearance of the second dust disk in the K band, where the first 
dust disk becomes transparent, would be easily blurred in real 
galaxies if perturbations from co-planarity occur, as has been 
observed in the form of "corrugations" in the Milky Way and 
external galaxies. 

We also investigated the approximation of the spiral arm 
component with an exponential disk (Sect. l4.2l l. and showed that 
this approximation has negligible effect on the predicted dust 
emission SEDs when the amount of dust is preserved. 

A particular feature of our model is the representation of the 
stellar emissivity SEDs in terms of just two parameters, S FR 
and old, since this enables UV/optical SEDs to be dereddened 
without recourse to population synthesis models, thereby avoid- 
ing bias due to the age/metallicity-opacity degeneracy. This fea- 



ture entailed the approximation of a fixed spectral shape for the 
SEDs of the old and young stellar populations. We quantified 
this approximation showing it has a negligible effect on the dust 
emission SEDs (Sect. 14.31 ). 

We have derived updated and improved solutions for the 
edge-on galaxy NGC 891 (Sect.O, one of the galaxies from the 
set used to constrain the geometry of the model. Since the lumi- 
nosity of the old stellar populations has been fixed to the values 
derived by Xilouris et al. (1999) from the optimisation of the op- 
tical/NIR images, old = 0.792 and B/D = 0.33, only 3 free pa- 
rameters were needed to fit the dust emission SED of NGC 891, 
namely r^, SFR, F. For the integrated SED of NGC 891 the 

best fit solution was found for = 3.5, SFR = 2.88M /yr 
and F = 0.41, which corresponds to a total dust luminosity 
of L'fi = 9.94 x 10 36 W, of which L d * f { = 6.89 x 10 36 W 

dust ' dusl 

is emitted in the diffuse medium (69%). The results indicate 
that the diffuse component dominates the emission longwards of 
60 yum and shortwards of 20 /mi, while the localised dust emis- 
sion within the star-forming complexes dominates at intermedi- 
ate wavelengths (20 - 60 /mi). It is only at intermediate wave- 
lengths (20 - 60 /im) where the localised dust emission within 
the star-forming complexes dominates. The total dust emission 
of NGC 89 1 is predominantly powered by the young stellar pop- 
ulations, which contribute 69% to the dust heating. The detailed 
input from the different stellar components is as follows: 11% 
from the bulge, 20% from the old stellar disk, 38% from the 
young stellar disk and 31% from the star forming complexes. 
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Fig. A.l. The extinction curve for the dust model used in this pa- 
per (solid line), which is a Milky Way dust model. With dashed- 
dotted line we plotted the observed mean extinction curve of our 
galaxy (Fitzpatrick 1999). Also plotted are the two components 
of extinction, the model absorption curve (dotted line) and the 
model scattering curve (dashed line). 
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Fig. A.l. The extinction curve for the dust model used in this pa- 
per (solid line) together with the contributions to the extinction 
from the different dust compositions used in the model: Si (dot- 
ted line), Gra (dashed line), PAH (dashed-three-dotted line). As 
in Fig. IA.ll the observed mean extinction curve of our galaxy is 
plotted with dashed-dotted line. 



Appendix A: Calculation of the extinction 
cross-section and scattering phase function 

A prerequisite for the calculation of radiation fields and infrared 
emission is knowledge of the extinction cross-section of the pop- 
ulation of grains, and, since we consider anisotropic scattering in 
the radiation transfer calculations, the scattering phase-function 
needs to be known as well. 

The absorption cross-section of grains of composition i = 
{Si,Gra,PAH°,PAH + }, C abs ,i, is obtained by integrating the ab- 
sorption efficiencies Q a bs,i over the grain size distribution n{a): 



C a bs,i(A) 



na n(a) Q abs j(a, A) da 



(A.l) 



where C a b s ,i is given in units of [cm 2 H _1 ], a mm is the minimum 
grain size and a max is the maximum grain size. 

Similarly, the scattering cross-section of grains of composi- 
tion ;, C sca j, is obtained by integrating the scattering efficiencies 
Qscaj over the grain size distribution n(a): 



C sea, i(.A) 



— I na 2 n 



(a) Q.ua,i(a, A) da 



(A.2) 



Then, by summing over the grain composition ; we obtain the 
total absorption and scattering cross-sections, C a ts and C sca : 



C„bs(A) - / C a b s ,i(A) 



Csca(A) - C Jefl , i(A) 



(A3) 
(A.4) 



The extinction cross-section C ext is the sum of the absorption 
and scattering cross-sections: 



Ce.xr(A) - C a b s (A) + C sca (A) 



(A.5) 



We note that the extinction cross section C ext is defined as per 
unit H. In some applications it is useful to define a cross section 



per unit dust mass, which we denote here C™ 7 : 
C ext (A) 



jT """ (4/3)^ fl 3 n ( a ) p g (Xa 



(A.6) 



where p g is the density of the grain material and C™ ( is in units 
of cm 2 /g. 

In turn, C"' xl is related to the extinction coefficient K ext , as 
used in the mathematical prescription of the dust distribution 
from Eqs.|4]and|6]using: 



K ext (A,R,z)=p d « s t(R,z)xC™ xt 



(A.7) 



where pd us t(R, z) is the dust mass density at position (R,z) in the 
galaxy in units of g cirT 3 and K exl ( A, R, z) is in units of cm -1 . 

Figs. IA. ll and lA.2l show the resulting extinction curve of the 
dust model adopted here, together with the absorption and scat- 
tering components (Fig. lA.U and the components given by the 
different grain composition (Fig. IA.21 . As expected, the figures 
confirm that the model extinction curve fits well the observed 
mean extinction curve of our galaxy. 

The averaged anisotropy of the scattering phase function g 
needed in the radiative transfer calculation is obtained in a simi- 
lar manner to Eqs. |AT] \A2\ |A3] andlA4l 

gi(A) = na 2 n(a)Q sca ,i{a,A)Q phase j{a,A)da (A.8) 



8(A) = 

i 

where Q P h ase , ; is the anisotropy efficiency. 



(A.9) 



Appendix B: The library of attenuations of stellar 
light for the diffuse stellar components 

The second set of simulated data needed to fit the panchromatic 
SEDs is the library of attenuations in the UV/optical/NIR as a 
function of and i. As mentioned in Sect. 12.41 the revision of 
the dust model required a recalculation of the database for the 
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1 -cos(i) 




0.0 
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0.4 0.6 

1 -cos(i) 



1 .0 



Fig. B.l. The attenuation-inclination relation for the disk (top), 
bulge (middle) and thin disk (bottom) in the B band. In each 
panel the solid curves represent the attenuations calculated with 
the dust model used in this paper, incorporating a mixture of 
silicate, graphite and PAH molecules. The dotted curves rep- 
resent the corresponding attenuations calculated with the dust 
model used in Paper III, incorporating only a mixture of silicates 
and graphites. In each panel the 7 different curves represent the 
attenuation-inclination relations for different r^: 0.1, 0.3, 0.5, 
1.0, 2.0, 4.0, 8.0, with the values increasing from bottom curves 
to top curves. 



attenuation of stellar light, as presented in Paper III. The over- 
all concept and characteristics of the calculations are the same 
as in Paper III, but a small change, mainly in the zero point of 
the calculations, was apparent due to the change in the relative 
contribution of absorption and scattering to the total extinction. 
Here we only give an example of a comparison between the 
attenuation-inclination curves obtained using the new and the 
old dust model (Fig. IB. II ). 

The attenuation-inclination relations for disks show a sys- 
tematic change with inclination when changing the dust model. 
Thus the attenuation for the low inclinations is decreased more 
than for the high inclinations, with a tendency for the curves 
to converge at the edge-on inclinations. This means that the 
shape of the attenuation-inclination is steepened for the present 
dust model. The curves for bulges show the biggest offset when 
changing the dust model, but in most of the cases there is no 
change in the shape of the curves. As one can see the curves run 
almost parallel, except perhaps for the lowest values of opacity. 
The smaller change is seen for the thin disk component, where 
neither the shape nor the zero point are changed significantly. 

We also did some tests to quantify the effect of the change in 
the dust model to the overall energy balance. By integrating the 
attenuation over all angles we obtained an estimate of the total 
energy absorbed by dust in a galaxy. This absorbed energy was 
found to be on average 10% smaller for the attenuations calcu- 
lated using the new dust model than for those from Paper III. 



Appendix C: Formulation of composite attenuation 
of stellar light 

In Sect. 5.1 of Paper III a generalised formula was given, show- 
ing how the composite attenuation (that is the overall attenua- 
tion of an arbitrary combination of luminosity components from 
stellar populations in the young stellar disk, the old stellar disk 
and the bulge) can be derived from the library of attenuation of 
stellar light from the diffuse component and the attenuation of 
the clumpy component. At any wavelength, the composite at- 
tenuation depends on the relative luminosity of the three stellar 
components, which we have described in this paper in terms of 
the parameters S FR, F, old and B/D, which we used to describe 
the dust emission. Here we re-write the generalised expression 
for the composite attenuation (Eq. 16 from Paper III) for the spe- 
cific parameterisation adopted in this paper. 

At a given wavelength A, the composite attenuation AniA in a 
galaxy is given by: 



Am A = -2.5 log — 



(C.l) 



where L9 and La are the intrinsic and the apparent luminosity 
densities. The quantities l}\ and La can be further expressed as a 
summation of the corresponding quantities for the disk, thin disk 
and bulge: 



jj) _ jj), disk _|_ ^0, tdisk , jQ, bulge 
A A A A 

L A = If* + Lf sk + L b A " lse 



(C.2) 
(C3) 



The apparent and intrinsic luminosity densities for the disk, 
thin disk and bulge are related as follows: 
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Am' 



disk 



Lf sk = LY isk lQ~ 



2 5 _ jO, disk j^disk 



Am' 



Idisk 



L tdisk = L Y disk (l-Ff A )lO 23" 



r 0, tdisk , 



jj), tdisk j^tdisk 
A A 



Am 



bulge 



Y bulge 



r 0, bulge 



10 2.5 



Y 0, bulge i bulge 
L A A A 



(C.4) 
(C5) 
(C.6) 



where Am" , Am', * and Am A 



bulge 



are the attenuation values ex- 



pressed in magnitudes for the diffuse component in the disk, thin 
disk and bulge and A d A ' sk , A' A lsk and A b A ulge are the corresponding 



attenuations expressed in linear form. Using Eqs. [131 [151 [191 and 
Eqs.|C4l|C3]E6]we can rewrite Eqs.lC2landlC3las: 



L° A = oldL°% t + 



$ FR jjoung 

IMsyr- 1 A ' mi ' 



I%Z + (B/D)oldL 



old 
A, unit 



(C7) 



U = oldL°"AT + 



old tdisk , S FR jjounv 

yr^ 1 

old a bulge 



^A,unit A 1 \M . yr^ ^^'Und 

(B/D) old L° A a lmjl A , 
By making the notation: 



e A (SFR,old) : 



C 17 D 1 T y °" ng 

SFR 1 L A.unit 



lM yr' old L°! d . 

wy A, unit 



(C8) 



(C9) 



the composite attenuation from Eq. lC.ll becomes: 

Am A = -2.5 log A A (CIO) 

where 



Af k + (1 - F f x ) e A (SFR,old)A tdisk + (B/D)A A ulge 



\ + e A {SFR,old) + B/D 



(C.ll) 



For the case that old - in the optical and in the UV Eq. IC. 101 

becomes: 



Am A = -2.5 log(l - F fx) + Am] 



tdisk 



(C12) 



The use of the calibration factor F ca \ in our procedure means 
that in practice the equations describing the attenuation due to 
the diffuse dust illuminated by the young stellar disk need to be 
rescaled to accommodate different values of F than those used in 
the calibration. For this we need to use the correction factor for 
the diffuse component corr^iF) as defined in Eq.|29j and rescale 
Eq. IC.5l in a similar way to the formulation of the radiation fields 
in Sect. 12321 

Lf sk = L° { tdisk ( 1 - F cal f A ) corr d (F) A' disk (C. 1 3) 

In this case Eq. lC.l ll becomes: 

A A = [Af k + (1 - F cal f A )corr d (F) e A (S FR,old) A' disk + 

(BID)A h A lse ]l[\ + e A (SFR,old) + B/D] (C.14) 
and Eq. lC.12l becomes: 

Am A = -2.5 log(l - F cal f A )corr d (F) + Amf sk 



(C.15) 



Eq. IC.14l together with Eqs. lC.lOl and lC.lll are the analog of the 
original expression for the composite attenuation from Eq. 16 
in Paper III. Eq. IC.15l together with Eq. lC.12l are the analog of 
Eq. 17 from Paper III. 



Appendix D: How to use the model to fit UV/optical 
to infrared/submm SEDs 

The six physical parameters determining our model prediction 
offheSEDare: 

- -4, SFR , F, old, B/D and Q 
whereby 

- r£, SFR, F, old, and B/D 

determine £T°ifL the model prediction for the dust luminos- 
ity density in the IR/submm as given by Eq. [43] of this paper, 
whereas 

- r£, SFR, F, old, B/D and i 

determine the attenuation in the measured UV/optical/NIR emis- 
sion (whereby the dependence on S FR and old is a weaker de- 
pendence due to the dependence of the composite attenuation on 
the relative amplitudes of the young and old stellar populations, 
as described in Appendix O. This attenuation is given in mag- 
nitude form, Am,,, by Eqs. IC. 10IC. 1 ll and lC14l In the following 
it is convenient to express it in the linear form: 



A A (r f B , S FR, F, old, B/D, i) = 10 ( 



-0.4xA»i|) 



(D.l) 



We note that two of the six physical parameters, SFR and 
old, are extrinsic (that is, the quantity scales with the amount of 
material in an object). This is a consequence of our model galaxy 
having a fixed size, expressed in terms of the fixed reference 
scale length of the old stellar population in B-band of hf s ^ = 

hf sk (B) = 5.67 kpc. 

The dust emission SED of the diffuse component of a galaxy 
with a valua3 for ftf u *differine from h d,sk , will be: 

^ — 1 s fc s.rel 



fdiff 



- A ,d'us t (SFR,old,r f B ,F,B/D) = ^x 



j diff, model ^ 
" 'A, dust 



ymodel 



node! 



(D.2) 



where 

SFR = SFR mode 'x( 2 



old = old modd - " 2 



(D.3) 

where S FR model is the star-formation rate of the model galaxy 
having the reference size h <.., old' nodel is the normalised lumi- 

& s, re j ' 

nosity of the old stellar disk population of the model galaxy hav- 
ing the reference size h sk „ S FR is the real star formation rate 

& s, re j ' 



14 i is the inclination of the disk, as used in Paper III 

15 Throughout this paper hf sk is taken as the true B-band scalelength 
of the stellar population. We note however that in any dusty galaxy hf sk 
will not actually be a directly observable quantity, even if the distance 
to the galaxy is known. As described and quantified in Paper IV the ap- 
parent size of a galaxian disk obtained from photometric fits will differ 
from the true size due to the stronger attenuation of light at the centre of 
galaxies compared to the outer regions. Due to the appearance of size 
(expressed by f) as a quadratic term in Eq.|D7I] determining the ampli- 
tude of the SEDs for a given value of 5 FR and old, this effect can ap- 
preciably influence the solutions obtained from fitting SEDs. Therefore, 
if using the measured angular sizes of galaxies as a constraint in fitting 
the dust emission SEDs the apparent sizes should be converted into true 
sizes using the correction factors tabulated in Tables 1-5 of Paper IV. 
Similarly, the B/D ratio should be converted from an observed ratio, 
as determined in photometric fits to data, into the intrinsic ratio, using 
correction factors tabulated in Paper III. In principle, i, as derived from 
measurements of axial ratios, should also be corrected. 
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of the galaxy that we want to model, old is the real normalised 
luminosity of the old stellar disk population of the galaxy that 
we want to model, and 



udisk 



j^disk 
s,ref 



(DA) 



Eq. lD.2l expresses the fact that radiation field energy density, and 
hence the colours of the dust emission, varies according to sur- 
face density of luminosity. In cases where a galaxy is unresolved 
t is unconstrained by the data, and becomes a further free param- 
eter of the model. Since we may also not know the distance D to 
the galaxy, it is convenient to express the dust emission SED of 
the diffuse component from Eq. ID.2l as a flux density by dividing 
throughout by AnD 2 , to obtain: 



S?? l Z at (.SFR,old,T f B ,F,B/D)=± 



,h disk f , 

\ s, rej / 



xL diff, model (S FRmo deI^ ^model^ ^ F B / D) 



(D.5) 



Correspondingly, if optical data is available, we can require 



that 



J L A> mriri, S FR, F, old, B/D, i, D) dA = 

J optical 



-J 



optical 



(Lf k + Lf sk +L b f se )dA 



(D.10) 



After manipulation of Eqs. Q~2] [13] Q3] and [19] from 
Sects. l2T2l l2~3?2l and l2331 and Eq.lD~7l this leads to: 

old - oM^ n ,SFR,F,B,D,i) - V™C1 + FUNC2) 



FUNC3 



where 
FUNCl 



= AnD 2 x ( 

J optical 



c obs 
A, star 



A^^SFR^^ld^/DJ) 



dA(DA2) 



FUNC2 = — 



SFR 

1M Q yr~ 



f 

J optical 



j young 1 j 
L A,unit dA 



(D.13) 



where 6 ga i is the half angle subtended at the Earth by the actual 
B-band scalelength of the galaxy: 



Vgal 



I disk 
n s 

D 



(D.6) 



We note that, provided the galaxy is sufficiently resolved for 6 ga i 
to be measured, S /{JL, as expressed by Eq. lD.5l depends on the 
value of t (via S FR and old -cf. Eq. ID. 31 ). but is independent of 
the distance D. 

Our model is constructed such that the total emitted lumi- 
nosity Lx, star of UV and optical light powering the dust emission 
can be directly constrained from available measured apparent 
UV/optical spatially integrated fluxes, Sf" t , corrected for at- 
tenuation (expressed here in linear form by Eq. ID. II ) as a func- 
tion of r^g, S FR, F, old, B/D, i and distance D: 



(t{, S FR, F, old, B/D,i, D) 



->A, star\T g 



AnD 2 



c obs 
J A, star 



A A {r{,SFR,F,old,BID,i) 



(D.7) 



where, as outlined in Sect. 5 of Paper III and in Appendix ICl of 
this paper the dependence on SFR, old and B/D is for the opti- 
cal range only. Depending on the exact range of wavelengths for 
which S^r' is available, L^star can be integrated over wave- 
length to obtain constraints on the parameters SFR = S FR" m 
and/or old - old co " as a function of and i. Specifically, if UV 
data is available we can require that 



FUNC3 = (1 + B/D) X L°Z, (D.14) 

Note that old "" is a weak function of SFR as well as r B , 
B/D and i due to the need to take into account the contribution 
of optical photons by the young stellar population, as described 
in Sect. 12370 

We are now in a position to determine the physical model pa- 
rameters from a combined set of measured UV/optical flux den- 
sities S° A bs slar , measured at wavelengths / i' obs < s,al \ and IR/submm 
flux densities of the pure dust component (corrected for con- 
tamination by direct stellar light at short infrared wavelengths) 
^°X d f measure d at wavelengths A s,dust . To do so, we min- 
imise the function 



s 

iobs.dust 



( ^ obs 
A, dust 



x2(SFR,old,T f B ,F,B/D) = 
■S^JSFR^ld^FB/D) 



@~iobs, dust 



(D.15) 



subject, if UV data is available, to the constraint SFR = 
SFR co "(t^, F, i, D) from Eq. lD.9l and. if optical data is available, 
old = old cm {T f B , S FR, F, B /D, i, D) from Eq. |D~TT1 cr^, dust are 
the 1 <t uncertainties in the measurements and S a, dust is given by 



1 A,du 



s ,(SFR,old,T f R ,F,B/D) 



jdiff 
A, dust 



(SFR,old,T B ,F,B/D) + 



L l A$st(SFR,F) 



AnD 2 



(D.16) 



f L A , star {r f B ,F,i,D)dA= f Lf sk dA (D.8) 
Juv Juv 

Combining Eq. \Ul from Sect. 12331 and Eqs. [DjJ and ES this 
leads to: 



SFR SFR co "( T f R ,F,i,D) . 

- = AnD 2 x 



1 



lM Q yr 



1M yr~ 



j young 
L unit, UV 



I 

Ju\ 



? obs 



A^Fi) 



dA 



(D.9) 



where S ff, and L l °°f are given by Eq. ID. 51 and Eq. [41] re- 
spectively. In the case that the distance, the optical structure 
and orientation parameters D, 6 ga i, (, B/D and i are known 
the optimization problem posed by Eq. ID. 151 is reduced from 
the parameter set (S FR,old,r B ,F,B /D,£) to the parameter set 

(S FR,old,r B ,F). With just 4 parameters, this might potentially 
be solvable purely considering the dust emission data, bearing 
in mind the orthogonal effect of these parameters on the am- 
plitude and colour of the dust emission SEDs, (as described in 
Sect. [6]l if at least 4 data points are available well sampling the 
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whole MIR/FIR/submm range. However, this will often be the 
exception rather than the rule, and in any case the fit is primar- 
ily and more robustly constrained through the optical and UV 
measurements. This is firstly because, if both UV and optical 
measurements are available the number of the primary search 
parameters would be reduced from four to just two - t b and F. 
Secondly, and perhaps more significantly, the model would no 
longer have any degree of freedom in terms of scaling parame- 
ters, due to the fact that, as noted above, SFR and old are the 
only extrinsic parameters in the full parameter set. 

Below we illustrate how to use the optical and UV con- 
straints by giving a possible processing path for a galaxy with 
a known (spectroscopically determined) redshift, for which inte- 
grated UV and optical photometry were available, and for which 
B/D, i and 8 ga i are known from optical imaging. The parameter 
set to be determined is thus r^, S FR, old, F: 

- step i): choose a trial value for each of t b and F 

- step ii): set SFR to SFR"'" and old to old"'" from Eqs.lD~9l 
and lD.lll for the trial values of t b and F. 

- step Hi): find SFR modd and old nwdel from Eq. Q33] substi- 
tuting for £ as defined in Eq. ID. 41 The value of hf sk used 
in Eq. ID. 41 should be derived from the optically measured 
value using the correction factors tabulated as a function of ; 
and t b in Tables 1 -5 of Paper IV. 



Appendix E: Tables describing the fixed parameters 
of our model 



Table E.l. The geometrical parameters of the model, as defined 
in Paper III and Sect. 12.11 All length parameters are normalised 
to the B-band scalelength of the disk of the standard modal 
galaxy, as derived for NGC 891: hf sk {E) = h iisk 



s,ref 



5670 pc 



/7f k (B) 

hf sk (V) 
hf sk (I) 
hf sk (J) 
hf sk (K) 

zf* 

/-disk 
"d 

^d 

^tdisk 

_tdisk 

/Jdisk 

"d 

-Rlisk 

^d 

Re 

b/a 



1.000 
0.966 
0.869 
0.776 
0.683 
0.074 
1.406 
0.048 
1.000 
0.016 
1.000 
0.016 
0.229 
0.6 



0.387 



step iv): using the trial values of v B and F, together 
with SFR nwde < and old model and £ from step Hi) find 
old m " del , t£, F, B/D). 



t diff, model fv pnmodel 
A, dust ( SFR 

step v): substitute L d ^J^™ odel from step iv) into Eq. ID. 51 to 



compute Sj ,f / ust . 



Compute L l ° c f ust {SFR C0 '\F) (note the use 
of the extrinsic value of SFR = S FR co " here). Substitute the 



values for S 



diff 



and L 



local 



, , , en,., , , determined at the wavelengths 

A, dust A, dust fe 

of the IR/submm observations in Eq. ID.15l to determine %2 
for the trial combination of and F. 

step vi): repeat steps i) to v), until the combination of t b 
and F that minimises %2 is found. The values of S FR and 
old from steps ii) to vi) found for the pair of t b and F that 
minimises %2 are then best fit parameters in SFR and old. 



We note that in the case of edge-on galaxies UV data should 
not be used to constrain SFR, since typically only a few percent 
of the total UV disk luminosity will be seen, and the solution 
can be subjected to stochastic variations, because the received 
photons may only be emitted by a small numbers of star-forming 
regions. In this case the SFR is better constrained from the FIR 
emission, as modelled for NGC 891 (Paper I and Sect. [3}. 

Having found the best fit parameter set, the following two 
further steps can be made: 



Table E.2. The intrinsic spectral luminosity densities of the old 
and young stellar populations, as defined in Sect. 12. 31 for old = 1 
and S FR = 1 M yr _I , respectively, and tabulated at the sam- 
pling wavelengths of the radiation fields. 



old=\ SFR =lM e yr- 



j young 



A 


W/Hz 


W/Hz 


912. 




0.344 x 10 21 


1350. 




0.905 x 10 21 


1500. 




0.844 x 10 21 


1650. 




0.863 x 10 21 


2000. 




0.908 x 10 21 


2200. 




0.926 x 10 21 


2500. 




0.843 x 10 21 


2800. 




0.910 x 10 21 


3650. 




1.842 x 10 21 


4430. 


4.771 x 10 21 


2.271 x 10 21 


5640. 


9.382 x 10 21 


3.837 x 10 21 


8090. 


19.54 x 10 21 


5.734 x 10 21 


12590. 


72.20 x 10 21 


0.931 x 10 21 


22000. 


64.97 x 10 21 


0.728 x 10 21 


50000. 


12.58 x 10 21 


0.141 x 10 21 



- step vii: deredden the UV/optical/NIR spectrum using the 
fitted values of tL, F, and the measured B/D and i 



- step viii: apply a population synthesis modelling fit to the 
dereddened UV/optical spectrum from vii to find the star 
formation history. 
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Table E.3. The parameters of the dust model considered in this 
paper (see Sect. I2.41 >. For parameter definition see Weingartner 
& Draine (2001). 



C 



parameter 


value 




b c (atoms/H) 


6.0 x 10-' 


xO.93 


b c ,i 


0.75 




b c ,i 


0.25 




Ooi (A) 


4.0 




"02 (A) 


20.0 




(T\ 


0.4 






0.55 




c s 


9.99 x 10- 


12 x 0.93 


a t , g (tan) 


0.0107 




a c>s (jxm) 


0.428 




a s 


-1.54 




Pg 


-0.165 




Pg(g/cm 3 ) 


2.24 




Si 






parameter 


value 




c s 


1.00 x 10- 


13 x 0.93 


a,, s (pm) 


0.164 




a c s (pan) 


0.1 




a s 


-2.21 




p s (g/cm 3 ) 


3.2 






0.3 





Table E.4. The wavelength dependence of the fraction of pho- 
tons escaping from the clumpy component into the diffuse 
medium as described in Sect. 12.21 for F = F ca i = 0.35. 



A 


(1-Fadfi) 


A 




912. 


0.427 


1350. 


0.484 


1500. 


0.527 


1650. 


0.545 


2000. 


0.628 


2200. 


0.676 


2500. 


0.739 


2800. 


0.794 


3650. 


0.892 


4430. 


0.932 


5640. 


0.962 


8090. 


0.984 


12590. 


0.991 


22000. 


0.994 


50000. 


0.999 



